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EMAP5:  A  3D  HYBRID  FEM/MOM  CODE 

Y.  Ji  and  T.  Hubing 

Department  of  Electrical  and  Computer  Engineering 
University  of  Missouri-Rolla 
Rolla,  MO  65401 


Abstract~EMAP5  is  a  numerical  software  package 
designed  to  model  electromagnetic  problems.  It 
employs  the  finite  element  method  (FEM)  to  analyze  a 
volume,  and  employs  the  method  of  moments  (MoM) 
to  analyze  the  region  exterior  to  the  FEM  volume.  The 
two  methods  are  coupled  by  enforcing  the  continuity  of 
the  tangential  fields  on  the  interface.  Dielectrics,  lossy 
dielectrics  and  metal  objects  can  be  modeled  within  the 
FEM  volume.  Metal  objects  can  also  be  located 
exterior  to,  or  on  the  surface  of  the  FEM  volume.  A 
simple  input  file  translator  (SIFT)  is  provided  with 
EMAP5  that  allows  users  to  define  basic  input 
configurations  on  a  rectangular  grid  without  a  graphical 
interface.  It  is  also  possible  to  use  commercial  mesh 
generators  to  analyze  more  complex  geometries  of 
arbitrary  shape.  EM  APS  can  model  three  types  of 
sources:  incident  plane  waves,  voltage  sources  on 
metal  patches  and  impressed  current  sources  in  the 
finite  element  region.  This  paper  describes  the 
implementation  of  EMAP5  and  provides  four  examples 
to  demonstrate  the  range  of  configurations  that  EMAP5 
is  capable  of  modeling  using  the  SIFT  translator. 

L  INTRODUCTION 

Methods  used  to  solve  Maxwell’s  equations  can 
often  be  classified  into  one  of  two  categories.  One 
category  is  integral  equation-based  methods.  Integral 
equation-based  methods  generally  solve  for  unknown 
fields  or  currents  on  the  surface  of  homogeneous 
volumes.  Only  the  surfaces  are  meshed;  therefore  the 
number  of  unknowns  is  relatively  small.  Integral 
equation-based  methods  are  well  suited  for  modeling 
homogeneous  structures  in  unbounded  geometries. 

The  other  category  is  partial  differential  equation 
(PDE)  -based  methods.  PDE-based  methods  require  a 
volume  mesh.  The  resulting  number  of  unknowns  is 
usually  much  higher  than  that  obtained  from  integral 
equation-based  methods,  but  the  reduced  complexity  of 
the  solution  method  can  more  than  compensate  for  this. 
PDE-based  methods  are  well  suited  for  modeling 
complex  inhomogeneous  structures,  particularly  those 
in  bounded  geometries  (e.g.  a  waveguide  or  a  metallic 
enclosure). 

The  finite  element  method  (FEM)  is  a  powerful 
PDE-based  method.  It  requires  the  entire  region  of 
interest  to  be  meshed.  However,  each  mesh  element 
can  have  its  own  electromagnetic  properties  (e,  |X,  a). 


Thus,  FEM  has  an  inherent  capacity  to  analyze 
inhomogeneous  structures.  In  addition,  FEM  does  not 
require  a  uniform  mesh.  Various  types  of  elements,  for 
example,  bricks,  tetrahedra  and  prisms,  can  be  used  to 
discretize  arbitrary  geometric  configurations. 
Therefore,  FEM  is  well  suited  for  modeling  complex, 
inhomogeneous  structures.  A  significant  disadvantage 
of  the  FEM  method  is  the  inability  to  model  unbounded 
geometries  without  using  absorbing  boundary 
elements.  Absorbing  boundary  elements  are  not 
effective  when  they  closely  conform  to  the  shape  of 
long  thin  objects  (e.g.  cables).  When  modeling  objects 
with  one  dimension  significantly  longer  than  the  others 
using  FEM,  it  is  generally  necessary  to  fill  the  empty 
space  around  these  objects  with  many  elements  that 
transition  from  the  object  of  interest  to  the  absorbing 
boundary.  This  can  significantly  inflate  the  overall  size 
of  the  problem. 

The  method  of  moments  is  not  actually  an 
electromagnetic  modeling  technique,  but  rather  it  is  a 
technique  for  solving  complex  linear  equations  [1,  2]. 
However,  electromagnetic  modeling  techniques  that 
apply  the  method  of  moments  to  the  solution  of  electric 
or  magnetic  field  integral  equations  are  generally 
referred  to  as  method-of-moment  (MoM)  techniques. 
MoM  techniques  can  be  used  to  model  metallic 
(particularly  wire)  structures  in  unbounded  geometries 
very  effectively. 

In  order  to  take  advantage  of  the  strengths  of  both 
the  FEM  and  MoM  methods,  several  researchers  have 
developed  hybrid  FEM/MoM  approaches  [3-6].  Hybrid 
FEM/MoM  methods  generally  take  advantage  of  the 
FEM’s  ability  to  model  complex  inhomogeneous 
structures  and  the  MoM’s  ability  to  model  unbounded 
geometries.  The  specific  strengths  of  a  particular  code 
depend  on  the  details  of  the  implementation, 
particularly  the  form  of  the  integral  equation  and  the 
basis  and  weighting  functions  employed. 

EMAP5  (ElectroMagnetic  Analysis  Program 
version  5)  is  a  hybrid  FEM/MoM  modeling  code 
developed  at  the  University  of  Missouri-Rolla  (UMR). 
EMAP5  was  developed  in  order  to  demonstrate  and 
validate  the  hybrid  technique  described  in  [6,7].  This 
technique,  used  in  conjunction  with  a  mesh  generator, 
is  particularly  useful  for  analyzing  EMI  from  printed 
circuit  board  geometries  with  enclosures  or  attached 
cables.  A  rudimentary  mesh  generator  (SIFTS)  is 
distributed  with  the  EMAP5  code.  This  non-graphical 


1054-4887  ©  2000  ACES 


2 


ACES  JOURNAL,  VOL  15,  NO.  1,  MARCH  2000 


user  interface  is  simple  and  portable,  making  it  a  good 
tool  for  educators  and  people  who  want  to  learn  to  use 
the  EMAP5  code.  However,  the  meshes  generated  by 
the  SIFT5  code  do  not  exploit  the  primary  advantage  of 
FEM  modeling  codes  (i.e.  the  ability  to  adapt  the  size 
and  shape  of  the  mesh  to  fit  the  geometry).  The 
EMAP5  and  SDFT5  source  codes  can  be  obtained  via 
the  web  at  <http://www.emclab.umr.edu/emap5>. 

H.  FORMULATION 


Although  details  of  the  EMAP5  formulation  are 
described  in  [6,7],  a  brief  summary  is  provided  below. 
The  general  structure  of  interest  is  shown  in  Figure  1. 
A  dielectric  volume  V2  has  electrical  properties  (e2,  (i2) 
and  is  enclosed  by  a  surface  S2.  A  conductive  volume 
V3  is  enclosed  by  a  conductive  surface  Sc.  The  fields 
within  V3  vanish.  Vj  denotes  the  volume  outside  of  V2 
and  V3,  which  is  assumed  to  be  free  space.  (Ei,  Hi)  and 
(E2,  H2)  denote  the  electric  and  magnetic  fields  in  Vi 
and  V2,  respectively.  The  unit  normal  vectors  for  S2 
and  Sc  are  defined  pointing  outward  toward  Vj.  The 
structure  is  excited  by  an  incident  wave  (Emc,  Hmc)  or 
impressed  sources  (Jmt,  Mmt).  The  scattered  electric  and 
magnetic  fields  are  (Es,  Hs).  The  objective  is  to  solve 
for  the  scattered  fields  (Es,  Hs)  or  the  surface  electric 
current  density  on  Sc. 


a 

n 


Figure  1.  A  dielectric  object  and  a  conductive  object 
illuminated  by  (E“c,  H“c )  or  (Jtat,  M“l). 


From  Maxwell  equations,  the  vector  Helmholtz 
equation  in  terms  of  E  can  be  written  as, 


Vx 


"VxE(r)" 


+  jco£0er  E(r) 


=  -  Jta,(r)  — 7 — - VxMtat(r). 

ja>t*o  Mr 


(1) 


After  multiplying  Equation  (1)  by  a  weighting  function 
w(r)  and  integrating  over  the  finite  element  domain  V2, 
one  obtains  the  FEM  weak  form  as  follows  [8,  9], 


Vll 


(VxE(r))«(7xw(r)) 


+  j  coe0£rE(r)  •  w(r) 


dV 


=  J  (nxH(r))«  vf(r)dS 


S  2 


jim(r)+ - VxMinl(r) 

JOMoMr 


w(r  )dV 


(2) 


Equation  (2)  can  also  be  derived  using  the  variational 
method  [10]  (pp.  158-161),  [11]  (pp.  158-160). 

Within  volume  V2,  EMAP5  uses  the  vector 
tetrahedral  elements  proposed  by  Barton  and  Cendes 
[10]  (pp.  56-59),  [12].  Each  basis  function  is  defined 
within  a  tetrahedron,  and  is  associated  with  one  of  the 
six  edges.  The  electric  field  E  within  volume  V2  can  be 
expanded  as, 

Nv 

E(r)  =  '£iEn  Wn(r)  (3) 

n=l 

where  Nv  is  the  total  number  of  interior  edges  and 
boundary  edges  on  the  dielectric  surface,  and  {En}  is  a 
set  of  unknown  complex  scalar  coefficients. 

The  surface  integral  term  in  Equation  (2)  can  be 
evaluated  by  using  a  surface  basis  function  f(r),  which 
is  discussed  later.  The  weighting  functions  chosen  are 
wn(r),  n=l...Nv.  After  multiplying  by  the  weighting 
functions,  Equation  (2)  can  be  discretized  as, 


'  An 

Aid 

'  Ei 

' 0 

o' 

"  O' 

\  inti 

=r 

+  [g  J 

_Adi 

Add. 

Ed. 

0 

Bdh. 

Jh. 

The  unknown  coefficients  [£„]  are  partitioned 
according  to  edge  type.  The  two  categories  are  interior 
edges,  which  are  denoted  by  the  subscript  i  in  Equation 
(4),  and  dielectric  boundary  edges,  which  are  denoted 
by  the  subscript  d  in  Equation  (4).  [ Jh ]  is  a  set  of 
unknown  complex  scalar  coefficients  for  the  surface 
electric  current  densities  on  Sd.  Sd  is  defined  as  S2  if 
the  conductive  body  is  not  adjacent  to  the  dielectric 
body;  otherwise,  Sd=  S2  -  (S2  n  Sc).  Edges  in  [Jh] 
include  dielectric  boundary  edges  and  junction  edges 
(i.e.  edges  that  are  on  both  the  conductor  and  dielectric 
boundary).  [Bdh\  is  a  square  matrix  if  there  are  no 
junction  edges  in  the  problem.  [gmt]  is  the  source  term, 
representing  sources  located  within  the  FEM  region. 
The  elements  of  [A],  [B&],  and  [gmt]  are  given  by, 
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T  (VxWn(r))»(VxWm(r)) 

Amn=  J  ja>MoMr  dV  (5) 

.  +  j  (0£0£r  w„(r)*  Wm(r)  . 

=  Jf.(r)-wB(r)dS  (6) 

Sd 

f  Jta+— - - (VxMta)  •wm(r)dV  (7) 

V2  [  JG>VoMr  J 

where  fn(r)  is  the  surface  basis  function.  Since  the 
tangential  E  field  vanishes  on  Sc,  Bmn  is  evaluated  over 
Sd  rather  than  all  of  S2.  A  detailed  description  of  how  to 
evaluate  the  terms  in  Equation  (5)  is  provided  in  [11] 
(pp.  254-257).  Appendix  A  describes  how  to  evaluate 
Equation  (6).  [gmt]  is  the  source  term  and  is  discussed 
in  a  later  section. 

The  MoM  surface  integral  equation  evaluated  by 
the  MoM  part  of  EMAP5  is  a  form  of  the  electric  field 
integral  equation  (EF3E)  [13], 

E“c  (r)  =  ^E(r)  +j  { M(r')  x  V'Go(r,r' ) 

+  j  koVoM^Goiry)  (8) 

-  j— V'«  J(r')  V'G0(r,r' 
k  o 


where  r  <=  S,  S  =  ScuS2,  7]o  and  k0  are  the  intrinsic 
impedance  and  wavenumber  in  free  space, 
respectively,  and 


Go(r,r')  = 


-ik  olr*”r1 


Ak  r-r'J 


(9) 


is  the  Green’s  function  in  free  space.  The  integral  term 
with  a  bar  in  Equation  (8)  denotes  a  principal-value 
integral.  The  singularity  at  r=r'  is  excluded.  The 
equivalent  surface  electric  and  magnetic  currents  are 
defined  as, 


J(r)  =  nxH(r),  reS 


GO) 


M(r)  =  E(r)  x  n ,  reSd.  (11) 

M(r)  vanishes  on  Sc.  J(r)  and  M(r)  can  be 
approximated  by  using  the  triangular  basis  function 
fn(r)  proposed  by  Glisson  [14].  On  the  surface  Sd,  the 
MoM  basis  function  fn(r)  and  the  FEM  basis  function 
wn(r)  are  related  by, 

Wn(r)  =  nxfn(r)  ,  r^Sa  (12) 


J(r),  M(r)  can  be  expanded  as, 


Ns 

J(r)=£/»f»(r)  (13) 

n=l 

Nd 

M(r)  =  ££»  f»(r)  (14) 

n=l 

where  Ns  is  the  total  number  of  edges  on  the  surface  S, 
and  Nd  is  the  total  number  of  edges  on  the  surface  Sd. 
{ Jn}  and  {En}  are  unknown  complex  scalar 
coefficients. 

The  weighting  functions  chosen  are  fn(r),  n=l,  ... 
Ns.  After  multiplying  by  the  weighting  functions, 
Equation  (8)  can  be  discretized  into  Equation  (15), 
which  is  a  matrix  equation.  Edges  on  Sd  and  Sc  are 
grouped  together  respectively. 


The  elements  in  matrices  [C],  [D]  and  [F1  ]  are  given 

by. 


Cmn=  -jkoTlo  Jfmfr)*  Jfn^Go^.r')^'  dS 

Sm  _Sn 


+  J— Jfm(r)*  -f(V'.fn(r'))V/Go(r,r')iS'  dS 


Dmn=  J  fm(r)< 


^  f  n(r')  x  V'G  o(r,  r')dS'+ — wn  (r)  dS 


Fm'=Jfm(r).Einc(r)dS.  (18) 

Sm 

Details  of  how  to  evaluate  Equation  (16)  are  provided 
in  [14].  Appendix  B  describes  how  to  evaluate 
Equation  (17).  Equation  (18)  can  be  evaluated  by  using 
the  Gaussian  Quadrature  technique. 

From  Equation  (15),  [7/,]  can  be  expressed  as, 


'o’]  [o  o  IF  o  1  r  o 

AIL0 

(19) 

where 

Cfi xh  —  Chh~~  ChcCcc  Cch 

(20) 

D'hd  -  Dm  -  Che  Cel  Dcd 

(21) 

K  =  ChcC:lFi-Fl  . 

(22) 
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After  combining  Equations  (4)  and  (19),  the  final 
matrix  equation  that  must  be  solved  is  given  by, 


"  W+ 

'0 

0 

\ 

~Ei 

\ 

0 

-BMC£D'hd_ 

y 

,Ei. 

=  [**]  + 


B+C?K 


(23) 


where  Xf  is  the  location  of  the  source  and  V  is  the 
voltage  magnitude,  and  S  (x)  is  the  Dirac  delta 
function.  Suppose  that  edge  m  coincides  with  the 
source.  The  normal  component  of  fm(r)  across  edge  m 
is  unity.  Thus,  the  source  term  given  by  Equation  (18) 
is, 

K=Vlm  (29) 


where  [A]  is  the  FEM  matrix  in  Equation  (4),  and 
BjhCfoD'hd  is  the  coupled  matrix  from  the  MoM 
equation.  The  right-hand  side  in  Equation  (23)  is  the 
source  term. 


The  far  fields  can  be  calculated  using  the 
equivalent  surface  currents  J  and  M  defined  in 
Equations  (10)  and  (11).  The  point  of  observation  is 
denoted  by  (r,  0,  cp)  in  spherical  coordinates.  Potential 
functions  A  and  F  can  be  approximated  by  [15], 


-;'*or 


A(r)  =  — - \j(r')ejk°T'msvdS’ 

Air  r  J  u 


4;rr 

-jlCQT 

e _ 

4;rr 


(24) 


-;*or 


F(r)  =  ^ - J  M(r')  eJkot'c05VdS' 


4ttt 

-;*o r " 


(25) 


4;rr 


where  r  cosy/ -  x  cos(p  sin#  +  y'sin^  sin#  +  z'cos#. 
EMAP5  evaluates  the  integral  terms  in  Equations  (24) 
and  (25)  using  the  Gaussian  Quadrature  technique.  The 
far  fields  are  given  by  [15], 

^  =  (-j'Mo)(4+— )  (26) 

^0 


where  lm  is  the  length  of  edge  m. 

Plane  wave  sources  are  also  described  using  the 
forcing  term  given  by  Equation  (18).  Assume  the 
incident  wave  can  be  expressed  in  a  spherical 
coordinate  system  as, 

E  ™=(E$e  +  E$)e-jk-r  (30) 

where  the  propagation  vector  k  is, 

k  =  £(sin  0  cos  (p  x  +  sin  6  sin  (p  y  +  cos#  z)  (31) 

and  (#,  (p)  defines  the  angle  of  the  propagation 
direction  of  the  plane  wave.  EMAP5  calculates  the  E 
field  of  the  plane  wave  on  the  surface  of  the  structure, 
then  evaluates  the  right-hand  side  of  Equation  (18) 
numerically. 

EMAP5  uses  current  filaments  on  tetrahedron 
edges  to  model  sources  located  within  the  FEM  region 
[11]  (pp.  324-325).  A  current  source  along  the  z-axis 
can  be  expressed  as, 

rl=lS(x-xf)S(y-yf)z  (32) 

where  (xf,  yf)  specifies  its  position,  I  denotes  the 
electric  current  magnitude,  and  d(x)  is  the  Dirac  delta 
function.  The  contribution  to  vector  [gmt  ]  in  Equation 
(7)  is  simply, 

g*"=ll\\  i»{w}8(x-xf)8(y-yf)dxdy  (33) 


E9=(-jW(\ (27) 
#0 

ffl.  SOURCE  AND  LOAD  MODELING 

EMAP5  can  model  voltage  sources  on  metal 
patches,  incident  plane  waves,  and  impressed  current 
sources  within  the  FEM  region.  A  delta-gap  voltage 
source  can  be  defined  anywhere  on  the  MoM  surface. 
It  can  be  modeled  as  an  incident  E-field  normal  to 
edges  coinciding  with  the  source  [9]  (pp.  260-261).  A 
delta-gap  voltage  source  parallel  to  the  x-axis  can  be 
expressed  as, 


where  l  is  the  edge  length  and  {W}  is  the  weighting 
function.  According  to  the  expression  in  Equation  (33), 
only  edges  coinciding  with  a  current  source  have  non¬ 
zero  components  in  the  [gmt]  vector.  After  the  E-fields 
along  the  current  source  edges  are  obtained,  the  voltage 
drop  across  the  source  can  be  calculated.  Thus,  the 
input  impedance  Zin  can  be  obtained.  The  desired  fields 
under  an  excitation  of  any  magnitude  voltage  Vin  can  be 
obtained  by  multiplying  the  factor  Vin/Zin  to  the  results 
obtained  when  1  is  set  to  one  Ampere. 

EMAP5  models  load  impedances  ZL  as  dielectric 
posts  on  tetrahedron  edges.  These  posts  have  a  finite 
conductivity  given  by, 


FT^VSix-x^x 


l 

ZLS 


(28) 


(34) 


Y.  Jl,  T.  HUB1NG:  EMAPS:  A  3D  HYBRID  FEM/MOM  CODE 


5 


where  l  is  its  length,  and  S  is  the  cross  sectional  area.  If 
the  load  is  treated  as  a  lumped  element,  its  contribution 
to  the  finite  element  matrix  is  as  follows  [11]  (pp.  324- 
325), 


The  new  solution  estimate  then  becomes, 

The  new  residue  and  bi-residue  are  calculated  as, 


k]=4-JJ  {W}{W}TS(x-xL)(y-yL)dxdy  (35)  rk+>  -  rk '  «kMp* 


where  (x^  yL)  is  the  position  of  the  impedance  load. 

IV.  MATRIX  EQUATION  SOLVER 

EMAPS  uses  the  complex  bi-conjugate  gradient 
method  (CBCG)  to  solve  the  final  complex  matrix 
equation,  Equation  (23).  The  major  advantage  of 
CBCG  lies  in  its  iterative  nature  and  its  ability  to  use 
the  row-indexed  method  to  store  sparse  matrices 
efficiently.  The  CBCG  method  was  first  introduced  by 
Jacobs  [16].  A  brief  outline  of  the  algorithm  is 
provided  in  the  following  section. 

A.  Complex  Bi-Conjugate  Gradient  Method 

Consider  a  complex  matrix  equation, 

[M][x]=[b]  (36) 

where  [M]  is  an  NxN  complex  matrix,  and  [b]  and  [x] 
are  order-N  complex  column  vectors.  Both  [M]  and  [b] 
are  known  while  [x]  is  unknown.  Four  column  vectors 
are  defined  as, 

1.  the  residue,  r, 

2.  the  bi-residue,  r, 

3.  the  direction  vector,  p, 

4.  the  bi-direction  vector,  p. 

These  vectors  are  initialized  as, 

ro=Po  =  b  (37) 

r0=Po=b*  (38> 

where  *  denotes  the  complex  conjugate.  The  Hermitian 
conjugate  of  M  is, 

Mh=(M*)t  (39) 

where  T  denotes  transpose.  The  inner  product  of  two 
complex  column  vectors  is  defined  as, 

<Pl.P2>  =(Pi*)TP2  .  (40) 

The  following  steps  are  repeated  N  times  at  most.  For 
k=0,  ...  ,  N-l,  calculate  the  step-length  parameter, 


a  =  <^5l> 

<  Pk>Mpk  > 


r*+i=rk-tf*MHP*  ’  (44) 

The  bi-conjugate  coefficient  is  calculated  using  the 
following  formula, 

R.  =.<bl!Pk’r^>  .  (45) 

<pk,Mp)t> 

The  new  direction  and  bi-direction  vectors  are  given 

by. 

Pk+i=rk+i+APk  (46> 

Pk+i=rk+i+A‘pk  •  (47) 

This  iteration  is  continued  until  a  termination  condition 
is  satisfied, 


B-<TOL 


where  TOL  is  a  number  defining  the  tolerance.  The 
best  value  of  TOL  depends  on  the  problem  for 
EMAP5.  Typically,  lxl0‘3  is  sufficient.  The  value  of 
TOL  is  defined  as  a  macro  in  the  EMAP5  source  code 
so  it  can  be  easily  adjusted  by  the  user. 

B.  Implementation 

It  is  not  necessary  to  explicitly  add  the  two 
matrices  on  the  left-hand  side  of  Equation  (23)  in 
CBCG.  The  FEM  matrix  [A]  is  a  sparse  and  symmetric 
matrix.  The  row-indexed  scheme  [17]  is  used  to  store 
it.  The  MoM  coupled  matrix  is  an  asymmetric 
compact  matrix.  A  fiill  2D  array  is  used  to  store  it. 
These  two  matrices  are  stored  separately  in  order  to 
make  it  is  easier  to  obtain  the  Hermitian  of  the  left- 
hand  matrix.  Since  [A]  is  symmetric, 

Ah  =  (AT)*  =  A*  (49) 

Thus,  it  is  not  necessary  to  search  the  non-zero 
elements  column  by  column  in  [A]  when  Hermitian 
operations  are  present.  Instead,  the  search  is  done  row 
by  row.  The  row-indexed  scheme  is  very  efficient  in 
this  case.  While  the  coupled  MoM  matrix  is  not 
symmetric,  exchanging  its  row  index  and  column  index 
will  result  in  its  transpose. 


(41) 


6 


ACES  JOURNAL,  VOL.  15,  NO.  1,  MARCH  2000 


El.sif 


El. far 


Pre-processor 


Field  Solver 


Post-processor 


Figure  2.  Components  of  the  EMAP5  package. 


C.  Accuracy  and  Convergence  Problem 

Although  the  CBCG  solver  conserves  memory,  it 
is  an  iterative  method  and  it  may  have  difficulty 
converging  in  some  situations.  Convergence  problems 
may  result  from  round-off  errors.  Usually  the  round¬ 
off  error  for  single-precision  data  is  1.1 9x1 O'7  [17].  If 
the  total  number  of  unknown  edges  is  1,000,  the 
average  error  due  to  round-off  errors  is 

Vl,000xL19xl0'7  =  3.8xlO'6. 

In  the  worst  case,  the  error  is  1, 000x1. 19xl0'7  = 
1.1 9x1  O'4.  Because  the  total  round-off  error  could  be 
very  large,  it  is  impractical  to  set  the  TOL  in  Equation 
(48)  to  be  on  the  order  of  1CT4  if  no  special  technique, 
such  as  pre-conditioning,  is  used  to  reduce  the  total 
round-off  error.  The  round-off  error  of  double¬ 
precision  data  is  2.22xl0‘16  [17].  If  double  precision  is 
used  for  the  matrix  calculations,  round-off  errors  will 
not  reduce  the  effectiveness  of  CBCG  considerably.  In 
EMAP5,  all  data  are  stored  in  double  precision; 
however,  double  precision  numbers  require  twice  as 
much  as  memory  to  store  as  single  precision. 

Another  potential  difficulty  with  the  CBCG 
solver  arises  when  the  algorithm  becomes  unstable. 
This  can  happen  when  the  denominator  of  the 
expressions  for  and  J3k  is  close  to  zero.  EMAP5 
monitors  the  value  of  the  denominator.  If  it  becomes 
too  small,  EMAP5  will  restart  the  algorithm  using  the 
current  solution  as  the  initial  estimate  as  suggested  by 
Jacobs  [16]. 

V.  COMPONENTS  OF  EMAP5 

EMAP5  is  written  in  ANSI  C.  Thus,  it  can  be 
compiled  without  modification  on  many  computer 
platforms.  As  shown  in  Figure  2,  the  EMAP5  software 
package  includes  three  major  components:  SIFTS, 
EMAP5  and  FAR.  EMAP5  was  developed  for  research 
and  educational  use.  It  does  not  have  a  sophisticated 
mesh  generator  or  graphical  visualization  tools.  SIFTS 
can  be  used  to  describe  simple  input  geometries  with  a 
rectangular  grid.  It  is  also  possible  to  use  commercial 


mesh  generators  to  analyze  more  complex  geometries 
of  arbitrary  shape.  EMAP5  is  the  FEM/MoM  field 
solver.  FAR  calculates  far-field  electric  fields  from  the 
equivalent  surface  currents  calculated  by  EM  APS. 

Standard  Input  File  Translator  version  5  (SIFTS) 
is  designed  to  generate  input  files  for  the  field  solver 
EMAP5.  SIFT5  reads  a  text  file  in  the  SIFT  format 
[18-20]  and  generates  an  input  file  for  EMAP5.  Users 
can  describe  the  structure  of  interest  by  using  eleven 
keywords.  The  physical  geometry  of  the  structure, 
source  type  and  location,  and  the  geometry  of  areas  in 
which  fields  are  to  be  printed  by  the  solver,  must  be 
specified. 

EMAP5  is  the  hybrid  FEM/MoM  field  solver.  It 
reads  an  input  file  generated  by  SIFTS.  A  .log  file  is 
generated  to  log  running  status  and  error  messages.  The 
log  file  was  initially  designed  to  debug  the  code,  but 
users  can  also  benefit  from  log  messages.  EMAP5  will 
print  out  fields  within  areas  specified  by  the  .sif  file,  to 
one  or  more  output  files.  All  equivalent  surface 
currents  J  and  M  must  be  printed  to  a  file  if  the  user 
needs  to  calculate  far  fields. 

FAR  is  a  program  used  to  calculate  the  far  field 
radiation  pattern.  FAR  needs  two  input  files.  One  is 
the  file  generated  by  SIFTS  and  the  other  contains  the 
equivalent  surface  currents  J  and  M.  The  FAR 
program  will  prompt  users  to  input  the  observing 
distance  R  in  wavelengths,  the  observing  6  interval  in 
degrees,  and  the  observing  <p  interval  in  degrees. 

VI.  STRUCTURE  OF  EMAP5 

Figure  3  shows  a  flowchart  for  the  field  solver 
EMAP5.  Descriptions  of  the  functions  used  by  EM  APS 
are  provided  below, 

•  ReadInputFile() 

reads  an  input  file  generated  by  SIFTS,  and 
checks  whether  the  input  file  is  consistent  (i.e.  the 
input  makes  sense).  If  not,  it  aborts  the  program 
and  reports  problems  with  the  input  to  the  user. 


Y.  Jl,  T.  HUBING;  EMAPS:  A  3D  HYBRID  FEM/MOM  CODE 


7 


•  ComputeBMatrix() 

places  a  3x3  sub-matrix  (corresponding  to  one 
observing  triangular  patch  and  one  source 
triangular  patch)  into  the  global  B  matrix  based 
on  the  local-to-global  mapping. 

•  ComputeCMatrixO 

places  a  3x3  sub-matrix  into  the  global  C  matrix 
based  on  the  local-to-global  mapping. 

•  ComputeDMatrixO 

places  a  3x3  sub-matrix  into  the  global  D  matrix 
based  on  the  local-to-global  mapping. 

•  ComputeCorrectionTerm() 

uses  B,  C,  and  D  matrices  to  compute  the  MoM 
coupling  matrix.  It  also  creates  the  necessary 
matrices  required  to  compute  the  equivalent 
surface  currents. 

•  FEMMatrixComputeO 

computes  a  6x6  matrix  for  each  tetrahedron  first, 
then  assembles  the  global  matrix  based  on  the 
local-to-global  mapping.  The  matrix  is  stored 
using  the  row-indexed  method. 

•  CreateRHSVector() 

creates  the  RHS  (right-hand  side)  vector  of  the 
final  matrix  equation. 

•  ConjugateSolver() 

uses  the  Complex  Bi-Conjugate  Gradient  Method 
to  solve  the  final  matrix  equation.  For  MoM-only 
analyses,  this  function  can  be  replaced  with  an 
LU  decomposition  solver. 

•  SurfaceFieldComputeO 

computes  the  equivalent  surface  currents  on  the 
boundary. 

•  PrintOutputO 

outputs  the  values  of  fields  in  regions  specified 
by  the  input  file  to  a  file  specified  by  the  input 
file.  All  equivalent  surface  electric  and  magnetic 
currents  must  be  written  to  a  file  if  users  are 
interested  in  the  far  field  pattern. 

VH.  NUMERICAL  RESULTS 

This  section  describes  four  example  geometries 
that  can  be  solved  using  the  EMAP5  software  package. 
Each  geometry  is  described  in  sufficient  detail  to  allow 
others  to  duplicate  the  analysis.  Results  obtained  from 
EMAP5  are  compared  to  results  obtained  using  other 
well-established  codes,  measurements,  or  analytical 
results. 


Figure  3.  EMAP5  flowchart. 


The  first  example  analyzes  the  electromagnetic 
scattering  from  a  square  plate  over  a  dielectric  slab.  As 
shown  in  Figure  4,  the  configuration  consists  of  a 
square  plate  of  dimension  0.75A  x  0.75A,,  and  a 
dielectric  slab  of  Er  =  2.0,  and  dimensions  of  0.4A,  x 
0.4A,  x  0.2A,.  The  plate  is  symmetrically  placed  above 
the  slab.  The  distance  between  the  plate  and  the 
dielectric  slab  is  0.1 5 A,.  The  system  is  illuminated  by  an 
x-polarized,  +z  traveling  plane  wave.  In  EMAPS’s 
hybrid  solution,  the  dielectric  is  divided  into  320 
tetrahedra.  The  surface  of  the  dielectric  and  the  metal 
plate  are  divided  into  392  triangle  patches.  Figure  5 
shows  the  normalized  far  field  pattern  calculated  by 
EMAP5.  For  comparison,  the  results  obtained  by  Arvas 
et  al.  [21]  have  been  plotted  on  the  same  scale.  The 
authors  in  [21]  used  a  method  of  moments  technique, 
which  is  based  on  a  combined  field  integral  equation 
formulation.  As  is  evident  from  the  figure,  the  two 
solutions  compare  very  well. 
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Figure  4.  A  conducting  plate  over  a  dielectric  slab. 


Figure  5.  The  normalized  far  field  pattern  for  the 
configuration  in  Figure  4. 


The  second  example  demonstrates  the  application 
of  the  EMAP5  code  in  modeling  microstrip  geometries. 
The  geometry  of  the  structure  is  shown  in  Figure  6. 
The  board  is  made  of  a  dielectric  with  e^.O.  The  trace 
is  excited  by  a  source  at  one  end,  and  is  terminated  by  a 
resistor  at  the  other  end.  Microstrip  lines  support  quasi- 
TEM  fields  [22]  (pp.  160-161)  and  can  therefore  be 
viewed  as  transmission  lines.  The  input  impedance  Zin 
of  a  transmission  line  is  given  by, 


(^L„  =  ;Z0tan(y3/).  (51) 


When  the  load  side  is  open,  the  input  impedance  is 
given  by, 


fein  \pen  J 


tan(yff/) 


(52) 


Thus,  the  characteristic  impedance  is  given  by, 
Z0=J(Zm)slw„{zJopm  .  (53) 

In  EMAP5’s  hybrid  solution,  the  dielectric 
volume  is  divided  into  900  tetrahedra.  The  surface  is 
divided  into  828  triangles.  Table  1  lists  (Zin)open, 

(Zin) short  ^  Z0  obtained  using  EMAP5  at  three 
different  frequencies.  The  results  in  Table  1  are 
consistent  and  Z0  =  56.4  Ohms  at  each  frequency. 
Figure  7  shows  the  input  impedance  at  the  source  end 
obtained  by  EMSIM  [23]  when  the  trace  is  terminated 
with  a  56.4  £2  resistor.  It  is  evident  that  the  microstrip 
line  is  almost  perfectly  matched  confirming  EMAPS’s 
calculation  of  the  characteristic  impedance. 


20  mm  4  mm  30  mm 


Figure  6.  A  microstrip  geometry,  (a)  y-z  plane  view, 
(b)  x-z  plane  view,  and  (c)  3D  view. 


7  —  7  J  Zq  0 

0  ZQ  +  jZLtan(0l)  (50) 

where  ZL  is  the  load  impedance,  (3  is  the  propagation 
constant;  /  is  the  length  of  the  transmission  line.  To 
determine  the  characteristic  impedance  Zo  of  the 
transmission  line,  one  can  find  the  input  impedance 
when  the  load  side  is  shorted  or  open,  respectively. 
When  the  load  side  is  shorted,  the  input  impedance  is 
given  by, 


Table  1.  Impedances  obtained  using  EMAP5. 


Source 

Frequency 

(MHz) 

( 7 .  ) 

\  in  /open 

(Ohms) 

(Zin  )short 

(Ohms) 

Zo 

(Ohms) 

300 

-j  1 87 

+jl7 

56.4 

500 

-j  1 06 

+j30 

56.4 

700 

-j69 

+j46 

56.3 
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Figure  7.  The  input  impedance  obatined  using 
EMSEM  with  a  56-Ohm  load  resistor. 


The  third  example  simulates  the  coupling 
between  an  active  trace  and  an  off-board  trace  on  a 
printed  circuit  board.  As  shown  in  Figure  8,  the  active 
trace  is  excited  by  a  1-volt,  600  MHz  source  at  one 
end,  and  is  terminated  by  a  50-Ohm  resistor  at  the  other 
end.  The  off-board  trace  runs  parallel  to  the  active  trace 
and  extends  well  off  the  board.  In  EMAPS ’s  hybrid 
solution,  the  dielectric  volume  is  divided  into  900 
tetrahedra.  The  surface  is  divided  into  876  triangles. 
For  the  first  test  case,  there  is  no  dielectric  material 
between  the  traces  and  the  ground.  Therefore,  the 
model  contains  only  perfect  electrical  conductor  (PEC) 
and  can  be  analyzed  using  only  the  MoM  portion  of 
EMAPS.  Figure  9  shows  the  current  on  the  off-board 
trace  obtained  by  EMAP5  compared  with  results 
obtained  by  EMSIM.  For  the  second  test  case,  the 
dielectric  constant  of  the  slab  is  set  to  4.0.  Figure  10 
shows  the  results  obtained  by  EMAP5  compared  with 
results  obtained  using  EMSIM.  For  both  cases,  good 
agreement  is  obtained  between  EMAP5  and  EMSIM. 

The  fourth  example  models  a  printed  circuit  board 
power  bus  structure.  As  shown  in  Figure  11,  the  top 
and  the  bottom  planes  of  the  board  are  PEC.  The 
dielectric  between  the  top  and  bottom  planes  has  a 
relative  permittivity  of  4.3(l-j0.02).  Two  test  ports  are 
identified  at  the  locations  shown  in  Figure  11.  In 
EMAPS ’s  hybrid  solution,  the  dielectric  volume  is 
divided  into  1,500  tetrahedra.  The  surface  of  the 
dielectric  and  the  metal  plates  are  divided  into  1,340 
triangle  patches.  The  model  is  verified  using 
measurements  of  an  actual  printed  circuit  board  with 
the  dimensions  shown  in  Figure  11.  Two  85-mil  semi- 
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Figure  8.  A  grounded  printed  circuit  board  with  active 
and  passive  traces,  (a)  y-z  plane  view,  (b)  x-z  plane 
view,  and  (c)  3D  view. 


Figure  9.  Current  distribution  on  the  off-board  trace 

(£r=1.0). 


Figure  10.  Current  distribution  on  the  off-board  trace 
te=4.0). 
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rigid  coaxial  probes  are  attached  at  the  port  locations. 
A  network  analyzer  is  used  to  measure  the  scattering 
parameters  of  the  two-port  system  up  to  2.0  GHz.  The 
measured  and  modeled  ISnl  and  IS2il  results  are  plotted 
in  Figure  12  and  Figure  13,  respectively. 

PEC 


jSi  i|  for  a  solid  power  bus  structure 


Figure  12.  The  measured  and  calculated  ISnl  for  the 
power  bus  structure. 


The  power  bus  structure  can  be  analyzed  as  a 
cavity  with  two  PEC  and  four  perfect  magnetic 
conductor  (PMC)  walls.  The  cutoff  frequencies  are 
given  as  follows  (for  lossless  case  and  |ir  =  1.0)  [22], 


mn 


m,n  —  0, 1, 2, ... 


(54) 


where  a  and  b  are  the  length  and  width  of  the  cavity, 
respectively;  m  and  n  are  the  mode  indices;  c  is  the 
light  speed  in  free  space;  er  is  the  relative  permittivity 
of  the  material  in  the  cavity.  For  this  power  bus 
structure,  only  TMZ  modes  are  excited. 


Table  2  lists  the  calculated  cutoff  frequencies  of 
the  first  five  TMZ  modes.  As  shown  in  Figure  12,  the 
ISnl  response  dips  at  these  resonant  frequencies.  At 
resonance,  more  power  is  delivered  to  the  board  hence 
more  power  is  transmitted  to  Port  2  as  shown  in 
Figure  13.  The  measurements,  numerical  results  and 
theoretical  analysis  agree  very  well. 


|S2i|  for  a  solid  power  bus  structure 


Figure  13.  The  measured  and  calculated  IS2il  for  the 
power  bus  structure. 


Table  2.  The  TMZ  mode  cutoff  frequencies  of  the 
power  bus  structure. 


Mode 

TMZ 

(1,0) 

TMZ 

(0,1) 

TMZ 

(1,1) 

TMZ 

(2,0) 

TMZ 

(2,1) 

fc 

(GHz) 

0.702 

1.019 

1.237 

1.405 

1.735 

Vm.  SUMMARY 

The  EMAP5  hybrid  FEM/MoM  code  allows 
users  to  analyze  geometries  with  large  PEC  surfaces 
and  lossy  dielectrics.  Unlike  codes  that  simply  use  the 
MoM  surface  as  an  absorbing  boundary  for  the  FEM 
region,  EMAP5  can  model  PEC  surfaces  inside, 
outside  and  adjacent  to  the  boundary  between  FEM  and 
MoM  regions.  This  permits  a  variety  of  interesting 
geometries  to  be  modeled  efficiently.  The  SIFTS  input 
translation  code  allows  users  to  input  simple  3D 
geometries  to  EMAP5  without  a  graphical  user 
interface.  Four  sample  geometries  were  analyzed  that 
demonstrate  the  types  of  problems  that  can  be  analyzed 
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using  SIFT5  and  EMAP5.  Good  agreement  was 
achieved  between  EMAP5  results  and  other  well- 
established  codes,  analytical  results,  and 
measurements. 
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Appendix  A.  Evaluation  of  B  Matrix 

The  elements  in  B  are  given  by. 


Jwm(r)*fn  (r)dS. 

(A-l) 

J±lm' 

Sm 

W  JI 

=  Jfm(r)«  ^^(rOxV'GoCr.r')^'  dS 

Sm  _Sn  _ 


When  edge  m  and  n  do  not  lie  on  a  common  triangle, 
Bm j  is  zero.  In  addition,  B ^  0  when  m=n,  because 

wn(r)  =  fixfn(r).  (A-2) 


If  edge  m  and  n  lie  on  a  common  triangle  Tq ,  Bmn  is 
given  by, 

Bmn  =  Jwm(r)*fn(r)<iS,=  J[nxfm(r)]  »fn(r )dS 

Sm  Sm 

=  Jn.[fm(r)xfD(r)]  dS  =  n.  J[fm(r)xfn(r)]  dS 

Sm  Sm 

=  Jk,  Xr»  +  <r«  -rm)^]dS  .  (A-3) 


The  ±  signs  before  lm  and  ln  refer  to  the  edge  directions 
of  edge  m  and  edge  n ,  respectively.  The  integral  term 
in  Equation  (A-3)  can  be  evaluated  using  the  one-point 
Gaussian  Quadrature  technique, 

B-  =  ^T^fi#[r”xr“+(r»-r«>)XrCP]  W 

4  AH 

where  rcp  is  the  centroid  of  triangle  Tq  . 


|(r'-r>V'Go(r,r')<iJ'  dS 


V'G0(r,r')  =  1  +  ^  ^°R  g~;1c°RR  (B-3) 

4?rR 

where  R  =  r  -  r',  R  =  |  R  | .  The  ±  signs  before  lm  and 

4  refer  to  the  edge  directions  of  edge  m  and  edge  n, 
respectively.  Because  the  EFDE  given  by  Equation  (8) 
already  excludes  the  singularity  at  r  =  r' ,  D'^  can  be 
evaluated  directly  using  the  seven-point  Gaussian 
Quadrature  technique. 


Appendix  B.  Evaluation  of  Matrix  D 

The  elements  in  matrix  D  are  given  by, 

J fm(r)»  -f fnCOxV'GoCr.rVS'  +  J-wn(r) 

5m  U  ^ 

=  -\Bm  +  ( since  B^-B^  (B-l) 

where  D'  is  defined  as  follows, 

m/i  7 
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Abstract 

A  hybrid  technique  is  developed  to  allow  the  scat¬ 
tering  from  small  appendages  to  be  approximately 
combined  with  the  scattering  from  a  large  body  of  rev¬ 
olution  (BOR).  The  hybrid  technique  thus  enables  the 
rotational  symmetry  of  the  large  BOR  to  be  exploited 
to  solve  the  scattering  problem  with  much  less  com¬ 
putational  complexity  than  a  fully  three-dimensional 
(3-D)  solution.  The  technique  combines  the  finite 
element  method  (FEM)  for  BOR  scattering  with  the 
method  of  moments  (MoM)  for  small  appendages. 
The  hybrid  formulation  is  discussed  in  detail,  includ¬ 
ing  an  approximate  Green’s  function  which  permits 
the  decoupling  of  the  solutions  from  the  FEM  and  the 
MoM.  Numerical  examples  are  given  to  show  the  ap¬ 
plicability  and  accuracy  of  the  hybrid  technique. 


1  Introduction 

The  symmetry  present  in  body  of  revolution  (BOR)  electro¬ 
magnetic  scattering  problems  permits  an  efficient  numerical 
solution  using  a  two-dimensional  (2-D)  technique  [1] — [51. 
However,  in  many  practical  problems,  the  rotational  symme¬ 
try  is  broken  by  the  presence  of  small  appendages  (see  Fig.  1). 
Thus,  a  three-dimensional  (3-D)  computational  method  is  re¬ 
quired  to  rigorously  compute  the  electromagnetic  scattering. 
Because  of  the  increased  computational  complexity  of  a  3-D 
method,  a  hybrid  method  is  developed  to  allow  the  scattering 
from  small  appendages  to  be  approximately  combined  with 
the  scattering  from  a  large  BOR  in  a  manner  similar  to  other 
hybrid  techniques  [6]-[9].  The  rotational  symmetry  can  then 
be  exploited  in  the  computation  of  the  scattering  from  the 
large  BOR.  The  hybrid  method  makes  use  of  the  finite  ele¬ 
ment  method  (FEM)  for  BOR  scattering  as  described  in  [4], 
and  the  method  of  moments  (MoM)  for  small  appendages  as 
described  in  [8].  The  remainder  of  this  paper  discusses  the 
hybridization  of  the  two  methods.  The  formulation  is  pre¬ 
sented  in  Section  2,  numerical  results  are  given  in  Section  3, 
and  concluding  remarks  are  found  in  Section  4. 


2  Formulation 

As  mentioned,  the  two  methods  to  be  hybridized  are  the  FEM 
for  BOR  scattering  and  the  MoM  for  small  appendages.  The 
FEM  for  BOR  scattering  makes  use  of  edge-based  vector  ba¬ 
sis  functions  to  expand  the  transverse  field  components  and 
node-based  scalar  basis  functions  to  expand  the  angular  com¬ 
ponent.  This  mixed  edge-node  formulation  eliminates  the 
problem  of  spurious  modes  [4],  and  it  directly  computes  all 
components  of  the  electric  field.  The  coupled  azimuth  po¬ 
tential  formulation,  which  computes  the  angular  component 
of  the  electric  and  magnetic  fields,  is  less  accurate  than  the 
mixed  edge-node  formulation  [5].  Mesh  truncation  for  the 
FEM  is  accomplished  using  cylindrical  perfectly  matched 
layer  (PML),  which  is  efficient  and  accurate.  PML  does  not 
alter  the  sparsity  of  the  FEM  matrix,  and  it  avoids  the  wasted 
computation  required  by  a  spherical  mesh  boundary  for  an 
elongated  geometry.  Finally,  the  FEM  equations  are  solved 
using  band  matrix  decomposition  techniques.  Validation  ex¬ 
amples  showing  the  efficiency  and  accuracy  of  the  method  are 
found  in  [4]. 

The  MoM  used  with  the  hybrid  method  is  straightfor¬ 
ward  and  makes  use  of  the  well-known  Rao-Wilton-Glisson 
(RWG)  basis  functions.  The  MoM  matrix  equations  encoun¬ 
tered  are  typically  of  low  dimension  and  are  thus  solved  using 
simple  LU-decomposition  techniques.  If  necessary,  fast  mul¬ 
tipole  techniques  can  be  employed  for  larger  problems  [10]. 

The  hybrid  formulation  is  best  understood  by  first  con¬ 
sidering  the  computation  of  the  scattering  from  a  BOR  with 
appendages  (see  Fig.  1)  entirely  by  the  MoM.  Assuming  that 
the  entire  target  consists  of  perfect  conductors,  the  applica¬ 
tion  of  the  MoM  uses  the  integral  equation 

E(r)  =  E*(r)  -  jk0m  j f  G0(r,r')  •  J(r ')dS'  (1) 

where  ko  =  is  the  free-space  wavenumber,  tjq  = 

y/fiojeo  is  the  impedance  of  free  space,  Go  is  the  free-space 
dyadic  Green’s  function,  E*  is  a  known  incident  electric  field, 
J  is  the  unknown  current  on  the  surface  of  the  target,  and  S 
is  the  surface  of  the  entire  target.  When  using  the  MoM,  J 


1054-4887  ©  2000  ACES 


14 


ACES  JOURNAL,  VOL  15,  NO.  1,  MARCH  2000 


Figure  1:  Example  of  a  large  BOR  with  small  appendages. 


is  found  by  discretizing  the  surface  S  and  applying  the  ap¬ 
propriate  boundary  condition  on  the  discretized  surface.  The 
scattered  field  is  then  found  from 

E5(r)  =  -jkon o  J ^  G0(r,  r')  •  J(r ')dS'.  (2) 

s 

This  3-D  method  generates  a  dense  matrix  equation;  there¬ 
fore,  if  the  geometry  is  large,  the  method  is  computationally 
very  intensive.  Further,  if  the  BOR  contains  inhomogeneous 
materials,  the  surface  integral  in  Eq.  (1)  must  be  replaced  with 
a  volume  integral,  further  increasing  the  computational  com¬ 
plexity. 

Because  of  the  unfavorable  computational  complexity  in 
solving  Eq.  (1)  with  the  MoM,  an  alternate,  hybrid  method 
is  sought.  In  the  hybrid  method,  the  MoM  is  applied  to  the 
small  appendages  only.  Thus,  Eq.  (1)  becomes 

E(r)  =  Ejj0R(r)  -  jk0r,0  J J  GBOr(t,t')  ■  3(r')dS'  (3) 

^App 

where  E|0R  represents  the  incident  field  on  the  appendages 
in  the  presence  of  the  BOR,  5app  represents  the  surface  of  the 
appendages,  and  Gbor  represents  the  dyadic  Green’s  function 
in  the  presence  of  the  BOR.  The  solution  of  Eq.  (3)  requires 
the  discretization  of  Sapp  rather  than  S  and  thus  is  much  less 
computationally  intensive.  The  incident  field  E ROR  is  calcu¬ 
lated  using  the  FEM  [4],  but  difficulty  is  still  encountered  in 
solving  Eq.  (3)  with  the  MoM  because  the  Green’s  function 
Gbor  is,  in  general,  unknown.  The  unknown  Green’s  func¬ 
tion  is  needed  not  only  for  the  solution  of  Eq.  (3),  but  also  to 
compute  the  scattered  electric  field  from 

Es (r)  =  -jkov o  J [  GB0R(r,r')  •  3(r')dS'.  (4) 

•S'app 

The  difficulty  posed  by  the  unknown  Green’s  function  Gbor 
is  alleviated  in  two  ways.  An  approximate  Green’s  function 
is  used  for  the  solution  of  Eq.  (3),  and  an  alternate  method 
based  on  the  reciprocity  theorem  is  used  to  compute  E5. 

There  are  four  steps  to  using  the  hybrid  method.  First,  the 
FEM  described  in  [4]  is  used  to  compute  scattering  from  the 
large  BOR  alone.  Then,  the  result  of  the  first  step  is  used 
to  compute  ER0R,  the  incident  field  on  the  appendages  in 
the  presence  of  the  BOR.  Next,  the  MoM  with  an  approxi¬ 
mate  Green’s  function  is  used  to  solve  Eq.  (3)  for  J,  the  cur¬ 
rent  on  the  appendages.  Note  that  in  the  usual  MoM  analy¬ 
sis  with  RWG  basis  functions,  the  unknowns  are  placed  on 


edges  which  are  shared  by  two  facets  in  the  mesh.  In  the  hy¬ 
brid  method,  unknowns  must  also  be  placed  on  the  boundary 
edges  connecting  the  appendages  to  the  large  BOR  to  account 
for  current  flow  from  the  appendages  to  the  large  BOR  [8]. 
Finally,  the  scattered  field  generated  by  J  is  evaluated  using 
reciprocity  and  added  to  the  scattered  field  from  the  BOR, 
which  is  found  in  the  first  step.  The  result  is  an  approxi¬ 
mation  to  the  scattered  field  from  the  entire  structure.  In  the 
remainder  of  this  section,  the  construction  of  the  approximate 
Green’s  function  and  the  computation  of  the  scattered  field  by 
reciprocity  are  each  discussed  in  turn. 


2.1  Approximate  Green’s  function 

The  approximate  Green’s  function  used  to  solve  Eq.  (3)  is  de¬ 
veloped  by  approximating  the  BOR  on  which  the  appendages 
reside  as  a  long  cylinder.  The  approximate  Green’s  function 
must  model  the  significant  field  interactions  between  points 
on  the  appendages  but  may  neglect  other,  less  significant  in¬ 
teractions.  The  interactions  are  classified  into  four  types.  The 
first  type  of  interaction  is  the  direct  path  interaction  between 
two  points,  and  examples  of  this  type  of  interaction  are  illus¬ 
trated  in  Fig.  2a.  The  second  type  of  interaction  involves  a 
reflection  by  the  large  BOR,  and  examples  are  illustrated  in 
Fig.  2b.  The  third  interaction  type  is  the  surface  wave  interac¬ 
tion  which  is  illustrated  in  Fig.  2c.  All  other  interactions  are 
classified  as  the  fourth  type.  Most  of  these  are  complex  inter¬ 
actions  which  are  not  modeled  by  the  cylinder  approximation 
to  the  BOR,  and  most  are  negligible. 

If  the  line  of  sight  between  two  points  on  the  appendages 
is  unobstructed  by  the  cylinder  approximation  to  the  large 
BOR,  the  first  two  interaction  types  can  be  computed  using 
the  half-space  Green’s  function.  This  is  done  by  finding  the 
reflection  point  on  the  cylinder  and  computing  the  tangent 
plane  to  the  cylinder  at  that  point.  If  the  line  of  sight  be¬ 
tween  the  two  appendage  points  is  obstructed  by  the  cylinder, 
both  of  the  first  two  interactions  are  zero.  For  many  prob¬ 
lems,  modeling  the  first  two  interactions  produces  acceptable 
results.  If  better  accuracy  is  desired,  the  third  type  of  interac¬ 
tion  can  be  computed  using  the  geometrical  theory  of  diffrac¬ 
tion  (GTD)  [11].  Also,  if  the  appendages  are  near  the  end 
of  the  BOR,  it  may  be  necessary  to  compute  an  interaction 
based  on  edge  diffraction  from  the  nearby  end.  Note  that  the 
effect  of  neglecting  some  of  the  interactions  is  to  neglect  the 
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corresponding  field  which  is  scattered  by  the  appendages,  re¬ 
flected  or  diffracted  back  to  the  appendages,  and  scattered  by 
the  appendages  again  [6]-[9]. 

In  all  of  the  examples  presented  in  this  paper,  the  approx¬ 
imate  Green’s  function  models  the  first  two  interaction  types 
(Fig.  2a-b)  and  neglects  all  other  interactions.  To  compute 
values  of  the  approximate  Green’s  function,  the  geometry  of 
the  protrusions  is  first  examined,  and  the  radius  of  the  infinite 
cylinder  used  to  approximate  the  large  BOR  is  found  from  the 
points  where  the  protrusions  attach  to  the  BOR.  Given  any 
two  points  on  the  protrusions,  the  value  of  the  Green’s  func¬ 
tion  is  then  computed  in  a  four  step  process.  First,  the  direct 
(line  of  sight)  path  between  the  two  points  is  computed.  If 
this  path  is  blocked  by  the  infinite  cylinder,  there  is  no  line  of 
sight  interaction.  Next,  the  interaction  involving  a  reflection 
by  the  infinite  cylinder  is  considered.  The  reflection  point  is 
found  using  Snell’s  law.  Note  that  it  can  be  shown  that  if  there 
is  no  direct  interaction,  there  is  also  no  reflected  interaction. 
Thus,  if  there  is  no  direct  path,  the  value  of  the  approximate 
Green’s  function  between  the  two  points  is  zero.  If  there  is 
a  direct  path,  the  third  step  is  the  application  of  a  coordinate 
translation  and  rotation  such  that  in  the  transformed  space,  the 
reflection  point  is  on  the  z  =  0  plane  and  the  normal  to  the 
cylinder  at  the  reflection  point  is  in  the  z  direction.  Finally, 
the  half-space  Green’s  function  is  applied  in  the  transformed 
space  to  compute  the  value  of  the  approximate  Green’s  func¬ 
tion  between  the  two  given  points. 

2.2  Computation  of  scattered  field 

While  the  use  of  the  approximate  Green’s  function  allows 
Eq.  (3)  to  be  solved  for  J,  the  approximation  is  not  accurate 
when  computing  E5  from  Eq.  (4).  Therefore,  an  alternate 
method  of  computing  Es  is  used.  The  alternate  method  is 
based  on  the  reciprocity  theorem.  Consider  a  short  dipole, 
located  at  point  r  and  oriented  in  the  u  direction.  From  the 
reciprocity  theorem, 

Es(r)  •  u  =  -jkoVo6 4/_r  J J  E£0R(r')  •  J(r')dS'  (5) 

Sapp 

where  E£0R  is  the  field  radiated  by  the  dipole  in  the  pres¬ 
ence  of  the  large  BOR.  Recall  that  this  field  can  be  computed 
by  the  FEM.  In  fact,  when  backscattering  is  being  computed, 
Ebor  same  as  Ebor  that:  is  used  in  Eq.  (3).  Thus,  all 
components  needed  to  compute  E5  using  Eq.  (5)  are  known. 

3  Numerical  Results 

Several  numerical  results  are  presented  to  show  the  validity 
and  capability  of  the  hybrid  technique.  In  all  of  the  results 
presented,  the  direct  and  reflected  interactions  are  modeled 
by  the  approximate  Green’s  function,  and  all  other  interac¬ 
tions  are  neglected. 


First,  the  validity  of  the  technique  is  tested  by  comput¬ 
ing  the  scattering  from  metallic  cylinders  with  capped  ends 
and  one,  two,  or  four  wings.  The  scattering  is  compared  with 
computations  from  the  Fast  Illinois  Solver  Code  (FISC)  [10], 
which  is  an  MoM  program  that  uses  a  multilevel  fast  multi¬ 
pole  algorithm  to  speed  up  the  matrix  solution.  The  cylinders 
considered  have  a  radius  of  1.25 A  and  height  of  5 A,  where  A 
is  the  ffee-space  electromagnetic  wavelength.  The  cap  at  one 
end  of  the  cylinder  is  pointed,  while  the  cap  at  the  other  end 
is  rounded.  The  attached  wings  are  1A  by  0.5A  by  0.025A. 


The  monostatic  RCS  as  a  function  of  azimuth  angle  from 
the  capped  cylinder  with  one  wing  is  shown  in  Fig.  3.  The 
agreement  between  FISC  and  the  hybrid  method  is  very  good. 
The  scattering  from  the  cylinder  without  the  wing  is  flat  as  a 
function  of  azimuth  angle,  so  it  can  be  seen  that  the  wing  has 
added  approximately  5  dB  peak-to-peak  swing  to  the  RCS 
in  the  VV-polarized  case  and  about  4  dB  peak-to-peak  swing 
in  the  HH-polarized  case.  Although  adding  a  second  wing 
presents  the  possibility  of  more  complex  interactions,  the  re¬ 
sults  shown  in  Fig.  4  continue  to  show  good  agreement  be¬ 
tween  FISC  and  the  hybrid  method.  The  peak-to-peak  varia¬ 
tion  in  the  RCS  is  about  4.5  dB  in  the  VV-polarized  case  and 
almost  6  dB  in  the  HH-polarized  case  when  two  wings  are 
present.  The  scattering  from  the  cylinder  with  four  wings  is 
shown  in  Fig.  5,  where  the  agreement  between  FISC  and  the 
hybrid  method  remains  very  good.  With  four  wings  present, 
the  peak-to-peak  swing  in  the  RCS  is  about  5  dB  in  the  VV- 
polarized  case  and  about  6  dB  in  the  HH-polarized  case. 


To  further  illustrate  the  capability  of  the  hybrid  method, 
two  more  computed  results  are  presented.  Both  of  these 
results  involve  1-GHz  scattering  from  a  missile  with  ap¬ 
pendages.  The  missile  is  12.5  m  (41. 7A)  long  and  has  a  radius 
of  0.625  m  (2.1A).  For  the  first  case,  a  3-cm  (0.1A)  by  3-cm 
(0.1A)  by  8.125-m  (27.7A)  ridge  is  located  on  the  missile  at 
azimuth  angle  0°,  and  computed  results  for  an  azimuth  scan 
and  for  an  elevation  scan  are  presented  in  Fig.  6.  Note  that  in 
the  azimuth  scan,  the  ridge  causes  the  scattering  to  vary  over  a 
3.5-dB  range  in  the  W-polarized  case  and  over  a  3-dB  range 
in  the  HH-polarized  case  while  the  missile  alone,  because  of 
its  rotational  symmetry,  has  a  constant  RCS  as  a  function  of 
azimuth  angle.  In  the  second  case,  two  fins  are  located  on  the 
missile  at  azimuth  angles  90°  and  -90°.  The  fins  are  trape¬ 
zoidal  in  shape  with  a  height  of  0.375  m  (1.25A),  bases  of 
1  m  (3.33A)  and  0.5  m  (1.67A),  and  a  thickness  of  0.01  m 
(0.03 A).  The  scattering  is  shown  both  for  an  azimuth  cut  and 
for  an  elevation  cut  in  Fig.  7.  In  the  azimuth  cut,  the  scatter¬ 
ing  has  changed  from  flat  for  the  missile  alone  to  a  function 
with  a  3  dB  peak-to-peak  variation  in  the  W-polarized  case 
and  a  2  dB  peak-to-peak  variation  in  the  HH-polarized  case. 
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'  Azimuth  Angle  (degrees)  Azimuth  Angle  (degrees) 

(a)  VV-pol  (b)  HH-pol 

Figure  3:  RCS  of  a  metallic  capped  cylinder  with  a  wing.  The  cylinder  has  a  radius  of  1.25A  and  height  of  5A,  and  the  wing  is 
1A  by  0.5A  by  0.025A. 


Azimuth  Angle  (degrees)  Azimuth  Angle  (degrees) 

(a)  VV-pol  (b)  HH-pol 

Figure  4:  RCS  of  a  metallic  capped  cylinder  with  two  wings.  The  cylinder  has  a  radius  of  1.25 A  and  height  of  5A,  and  the 
wings  are  1A  by  0.5A  by  0.025A. 


Azimuth  Angle  (degrees)  Azimuth  Angle  (degrees) 

(a)  W-pol  '  (b)  HH-pol 

Figure  5:  RCS  of  a  metallic  capped  cylinder  with  four  wings.  The  cylinder  has  a  radius  of  1.25A  and  height  of  5A,  and  the 
wings  are  1A  by  0.5A  by  0.025A. 
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Azimuth  Angle  (degrees)  Elevation  Angle  (degrees) 

(a)  Azimuth  cut  (b)  Elevation  cut 

Figure  6:  RCS  of  a  missile  with  a  ridge  at  1  GHz.  The  missile  is  12.5  m  (41. 7A)  long  and  has  a  radius  of  0.625  m  (2.1A).  The 
ridge  is  3  cm  (0.1A)  by  3  cm  (0.1A)  by  8.125  m  (27.7A)  and  located  at  azimuth  angle  0°. 


(a)  Azimuth  cut 


(b)  Elevation  cut 


Figure  7:  RCS  of  a  missile  with  two  fins  at  1  GHz.  The  missile  is  12.5  m  (41. 7A)  long  and  has  a  radius  of  0.625  m  (2.1A). 
The  fins  are  trapezoidal  with  a  height  of  0.375  m  (1.25A),  bases  of  1  m  (3.33A)  and  0.5  m  (1.67A),  and  a  thickness  of  0.01  m 
(0.03A);  they  are  located  at  azimuth  angles  90°  and  —90°. 
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4  Summary 

The  hybrid  method  is  a  useful  extension  of  the  FEM/BOR 
capability  described  in  [4].  While  rigorous  computation  of 
scattering  from  a  BOR  with  appendages  requires  a  3-D  com¬ 
putational  method,  the  hybrid  method  separates  the  BOR  part 
of  the  problem  from  the  appendages.  The  rotational  sym¬ 
metry  of  the  BOR  part  of  the  problem  can  then  be  exploited 
for  computational  efficiency  while  only  the  appendage  part, 
which  is  typically  much  smaller  than  the  BOR  part,  requires 
a  3-D  method.  Numerical  results  show  the  impact  of  the  ap¬ 
pendages  on  the  scattering  from  the  entire  structure  and  verify 
the  validity  and  capability  of  the  hybrid  method. 
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ABSTRACT 

A  study  has  been  carried  out  on  the  impedance  of  a  half¬ 
wave  dipole  above  a  finite  ground  plane.  The  Method  of 
Moment  code,  NEC-4,  has  been  used  and  it  has  been 
found  that  for  ground  planes  which  are  greater  than  1.0  by 
1.0  wavelengths,  the  changes  from  the  effect  with  an  infi¬ 
nite  plane  are  negligible.  With  very  small  planes,  the 
current  distribution  across  the  ground  plane  and  the 
impedance  are  substantially  different  from  those  with 
larger  ground  planes. 

1.  INTRODUCTION 

There  is  a  dearth  of  information  in  the  open  literature  on 
the  effect  of  a  finite  ground  plane  on  the  impedance  of  a 
nearby  halfwave  dipole  although  several  authors  have 
considered  the  effect  of  finite  sheets  on  radiation  patterns. 
We  examine  the  general  problem  of  impedance  as  a  func¬ 
tion  of  sheet  size  and  dipole  position  relative  to  the  edges 
of  the  sheet.  Dipole  impedances  have  been  calculated 
using  the  Method  of  Moments  Code,  NEC-4  [1],  for  var¬ 
ious  dipole  positions  over  square  sheets  of  side  0.5,  1.0, 
2.0  and  2.5X  and  are  examined  along  with  theoretical 
considerations. 


2.  COMPUTATIONS 

The  dipole  was  of  diameter  0.002  X  (2mm)  and  was  0.5X 
(500mm)  long.  Spacing  from  the  ground  plane  was 
0.221  A,  (221  mm).  Throughout  this  paper,  all  dimen¬ 
sions  are  quoted  in  mm  and  in  wavelengths  at  300  MHz. 
The  RF  performance  was  examined  over  the  frequency 
band  270-330  MHz. 

Square  sheets  of  size  0.5, 1.0,  2.0  and  2.5A,  (500  to  2500 
mm)  were  tested.  The  dipole  position,  initially  in  the  cen¬ 
tre  of  the  sheet,  was  altered  in  steps  of  0.25X,  (250  mm) 
in  either  the  x  or  y  directions.  The  dipole  lay  along  the  x 
axis  (Figure  1). 


Ay 


X 


"  Dipole 


Figure  1  Geometry  of  Dipole  over  a  Finite  Square 
Ground  Plane.  S  was  varied  between  0.5  and  2.5  X. 

The  ground  planes  were  modelled  as  a  square  mesh  with 
segments  which  were  0.05A,  long  (50  mm)  and  had  a 
radius  of  008  X  (8.0  mm)  giving  a  length  to  diameter  ratio 
of  3.125  which  is  greater  than  the  required  lower  limit  of 
2.0.  This  diameter  is  derived  from  the  equal  area  rule  but 
a  convergence  test  was  also  applied  (Section  7). 

The  dipole  was  modelled  as  21  perfectly  conducting  lin¬ 
ear  segments  in  a  length  of  0.5  wavelengths  (500  mm); 
each  segment  had  a  diameter  of  0.002  X  (2  mm).  The 
excitation  was  a  voltage  source  of  1  volt  applied  to  the 
centre  segment.  Trials  were  carried  out  with  51  segments 
for  the  dipole  in  free  space  but  the  results  were  very  sim¬ 
ilar  to  those  with  21  segments  (Section  7). 

Computations  were  carried  out  with  both  double  and  sin¬ 
gle  precision  versions  of  NEC-4  (Section  7). 


3.  FREE  SPACE 

An  initial  run  was  made  with  the  dipole  with  no  ground 
plane  present.  Table  1  gives  the  impedance  values,  the 
resonance  frequency  being  just  below  285  MHz  where 
the  length  would  have  been  0.473A,. 
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Belrose  [2]  NEC-4 


Figure  2  Reactance  against  Resistance  for  a  halfwave 
dipole  in  free  space  The  NEC -4  values  are  indicated  at 
5*0  MHz  steps  between  270  and  330  MHz  (Table  1). 


Table  1  Impedance  of  dipole  in  free  space  (Dipole 
0.002A,  diameter,  0*5X  long  at  300  MHz) 


F 

MHz 

NEC-4 

Belrose  [2] 

R 

(ohms) 

X 

(ohms) 

R 

(ohms) 

X 

(ohms) 

270 

61.1 

-46.9 

66.1 

-25.0 

275 

64.5 

-30.9 

280 

68.2 

-14.9 

285 

72.1 

1.02 

290 

76.1 

16.9 

295 

80.5 

32.8 

300 

85.0 

48.7 

87.6 

65.0 

305 

89.8 

64.6 

310 

94.9 

80.5 

315 

100.3 

96.5 

320 

106.0 

112.6 

325 

112.1 

128.7 

330 

118.5 

145.0 

117.4 

155.8 

Figure  2  compares  the  computed  impedances  and  those 
calculated  from  Equations  15.16  and  15.17  of  Belrose, 

[2]  which  are  for  a  monopole.  These  values  were  doubled 
to  obtain  the  numbers  for  a  dipole  shown  in  Figure  2  and 
Table  1.  The  agreement  is  good.  It  implies  that  the 
antenna  is  about  3%  longer  electrically  than  its  dimen¬ 
sions  would  suggest  from  Belrose.  In  their  Figure  13.31, 
showing  calculated  values  of  resistance  at  resonance  as  a 
function  of  wavelength/diameter,  Schelkunoff  and  Friis, 

[3] ,  give  a  value  of  64  ohms  for  a  dipole  of  the  dimen¬ 
sions  used  here  whereas  NEC-4  and  Belrose  both  give  72 
ohms. 

4.  GROUND  PLANE 

Figure  3  to  Figure  5  show  the  impedance  for  the  central 
and  the  two  extreme  positions.  The  values  are  indicated 
at  5,0  MHz  steps  between  270  and  330  MHz. 

For  sheets  of  side  1.0 A,  and  larger,  there  is  no  significant 
change  in  impedance  with  sheet  size  when  the  dipole  is 
centrally  placed  or  when  it  is  over  a  parallel  edge.  When 
the  dipole  is  moved  longitudinally  so  that  its  centre  is  on 
the  edge,  there  is  more  variation  with  sheet  size  as  Figure 
5  demonstrates.  Results  for  the  X/2  square  sheet  also 
show  significantly  lower  resistance  values.  The  imped¬ 
ance  values  for  an  infinite  ground  plane  have  been 
included  in  Figure  3  and  show  that,  for  plates  of  side 
greater  than  1.0  X,  the  results  have  converged.  For  the 
larger  sheets,  there  does  appear  to  be  a  very  slight  cyclical 
variation  with  dipole  position.  This  is  shown  in  Table  2 
for  300  MHz  on  a  2.5A-  square  sheet. 

- 0.5  1.0  2.0  -*•-  2.5  Infinite 


Figure  3  Impedance  of  dipole  (Central  position) 
between  270  and  330  MHz 
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0.5  1.0  2.0  2.5 


Figure  4  Impedance  of  dipole  (over  parallel  edge) 
between  270  and  330  MHz 


“•*-  0.5  1.0  ”J*‘"  2.0  ••'•25 


Figure  5  Impedance  of  dipole  (over  perpendicular 
edge)  between  270  and  330  MHz. 


Table  2  Impedance  of  dipole  above  ground  plane  of 
side  2.5A  as  a  function  of  position  off-centre.  (Fre¬ 
quency  300  MHz) 


y/X 

II 

o 

x/A 

y/X  =  0 

0 

93.9  +  j  90.2 

0 

93.9  +  j  90.2 

0.25 

94.0  +  j  90.2 

0.25 

94.1  +  j  90.3 

0.5 

93.9  +  j  90.1 

0.5 

93.8  +j  90.0 

0.75 

94.1 +j  90.1 

0.75 

94.7  +  j  90.5 

1.0 

92.6  +j  90.6 

1.0 

87.7  +  j  91.6 

1.25 

94.5 +j  77.8 

1.25 

77.8 +  j  66.0 

5.  CURRENT  FLOW 

As  noted  above,  the  values  of  the  impedance  for  the 
smallest  ground  plane  (0.5  by  0.5  A,)  were  quite  different 
from  those  of  the  larger  ground  plane. 

The  currents  on  the  ground  planes  at  intervals  of  0.25 
wavelengths  have  been  plotted  for  a  dipole  centrally 
placed.  The  abscissa  (current  magnitude)  has  been  multi¬ 
plied  by  10^  for  convenience. 


'  0.0  10.0  20.0  30.0  40.0  50.0 


Segment  No 


Figure  6  Current  Magnitude  (x  10^)  on  a  Ground 
Plane  parallel  to  the  dipole  and  immediately  below  the 
dipole 

Figure  6  shows  the  currents  on  the  ground  plane  immedi¬ 
ately  below  the  dipole  for  several  ground  plane  sizes 
while  Figure  7  and  Figure  8  show  the  currents  0.25  and 
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0.5  wavelength  away  from  the  centreline  underneath  the 
dipole.  The  most  notable  feature  is  that  the  current  distri¬ 
bution  is  remarkably  alike  for  ground  planes  with  sides 
1.0  to  2.5  wavelengths  while  the  ground  plane  of  side  0.5 
wavelength  is  very  different. 


- S-0.5  - S-J.O  - $-2.0  - S-2.5 


Figure  7  Current  Magnitude  (x  10^)  on  a  Ground 
Plane  parallel  to  the  dipole  and  offset  by  0.25  wave¬ 
length  from  the  dipole 


- S-I.O - S-2.0  —  S-2.5 


Segment  No 

Figure  8  Current  Magnitude  (x  10^)  on  a  Ground 
Plane  parallel  to  the  dipole  and  offset  by  0.5  wave¬ 
lengths  from  the  dipole 


6.  CHECKS 

It  was  possible  to  provide  some  corroboration  for  the  pre¬ 
dictions  by  computing  the  mutual  impedance  provided 
when  a  half-wave  dipole  is  over  an  infinite  ground  plane. 
The  mutual  impedance  is  then  the  impedance  of  the 
dipole  alone  minus  the  mutual  impedance  of  the  dipole 
over  an  infinite  ground  plane.  These  are  compared  with 
results  from  Smith  [5]  and  STC  Data  [4]  in  Table  3. 
These  results  are  in  better  agreement  with  the  STC  data. 

Since  the  image  of  the  dipole  in  the  ground  plane  is  neg¬ 
ative,  the  computed  figures  are  in  opposite  senses  to  those 
from  two  parallel  dipoles  at  the  same  spacing  as  the 
dipole  and  its  image.  The  computed  figures  have  been 
reversed  in  sign  for  comparison  with  other  sets  of  data. 


Table  3  Mutual  Impedance  of  a  Dipole  above  an  Infinite 
Ground  Plane 


Af 

MHz 

NEC-4 

STC  data 

R  A  Smith 

R 

jx 

R 

jx 

R 

jx 

-5.0 

1.1 

-31.4 

1.1 

-34.8 

4.0 

-36.9 

-1.0 

-1.5 

-33.6 

-1.5 

-36.4 

1.5 

-36.3 

+2.0 

-3.7 

-35.2 

-3.2 

-35.5 

-0.3 

-36.0 

+4.0 

-5.3 

-36.2 

-4.2 

-35.1 

-1.3 

-35.7 

+6.0 

-7.0 

-37.2 

-5.2 

-34.5 

-2.5 

-35.3 

+7.0 

-7.9 

-37.8 

-6.0 

-34.1 

-3.1 

-35.1 

RMS  Error  in  NEC 

1.2 

2.6 

3.8 

;  3.0 

7.  CONVERGENCE  TESTS 

Several  convergence  tests  were  carried  out  on  the  input 

data. 

1)  Single  and  double  precision  versions  of  NEC-4 
were  run  on  all  the  examples  for  a  ground  plane 
with  a  side  of  0.5  X.  The  changes  in  all  values  of 
impedance  were  less  than  0.001%. 

2)  The  radius  of  segments  in  the  ground  plane  was 
varied  between  20  and  2.5  mm  with  the  segment 
length  held  constant  at  50  mm  (Figure  9).  The 
value  chosen,  8.0  mm,  was  based  on  the  equal  area 
rule  and  computed  results  agreed  with  other 
results  to  within  1  or  2  ohms. 


3)  The  number  of  segments  in  the  dipole  was  varied 
from  21  to  51.  Differences  of  less  than  1.0  degree 
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in  phase  and  1%  in  impedance  magnitude  were 
observed. 


- 270.0  MHz - 300.0  MHz  - 330.0  MHz 


Figure  9  Impedance  for  three  frequencies  for  dipole  in 
centre  of  0.5  A,  ground  plane  with  the  segment  radii  of 
the  ground  plane  varied.  The  impedance  magnitude  is 
shown  normalised  to  the  value  at  a  segment  radius  of 
2.5  mm 

8.  CONCLUSIONS 

1)  The  impedance  of  an  isolated  0.5A  dipole 
obtained  by  NEC  is  closely  modelled  by  the  Bel- 
rose  [2]  formula  for  a  dipole  of  the  same 
dimensions.  The  difference  of  3%  in  frequency 
may  arise  from  small  errors  in  Belrose’s  determi¬ 
nation  of  electrical  length  which  was  empirical. 
Alternatively  it  may  represent  a  genuine  differ¬ 
ence  between  the  two  models.  It  has  been  noted 
(Section  7)  that  small  changes  in  the  NEC-4 
model  resulted  in  changes  in  the  impedance  of  up 
to  1%  in  magnitude. 

2)  When  mounted  above  the  centre  of  a  finite  ground 
plane,  the  dipole  impedance  does  not  alter  signif¬ 
icantly  for  sheets  of  lX  side  or  larger.  The 
resistance  of  a  dipole  above  a  sheet  of  side  0.5A  is 
significantly  lower.  The  plots  of  current  flow  in 
Figure  6  to  Figure  8  support  the  impedance 
results. 

3)  When  the  dipole  is  over  the  parallel  edge  of  the 
sheet,  the  general  picture  is  similar  to  that  of  the 
central  position. 

4)  When  the  dipole  centre  is  over  the  perpendicular 


edge,  there  is  more  variation  with  sheet  size. 

5)  The  impedance  on  larger  sheets  appears  to  be 
more  sensitive  to  movement  along  the  dipole  axis 
than  to  movement  in  the  perpendicular  direction. 
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Abstract 

Design  procedures  for  planar  waveguide  slot  arrays  have  been 
known  and  used  for  years,  but  the  inexperienced  designer  is 
often  faced  with  a  number  of  practical  problems  when  trying 
to  implement  them,  especially  in  the  case  of  larger  arrays 
consisting  of  sub-arrays.  Some  helpful  hints  concerning  the 
numerical  implementation  are  provided.  These  include  an 
effective  procedure  to  avoid  the  necessity  of  good  initial 
guesses  for  the  unknown  dimensions  and  recommendations 
for  the  subdivision  of  the  nonlinear  equations  into  smaller 
groups.  Other  practical  aspects  which  are  not  explicitly 
defined  elsewhere  in  the  literature,  are  also  addressed. 

1.  INTRODUCTION 

Planar  slot  arrays  are  key  components  of  high  performance 
airborne  radar  systems,  and  also  find  application  in  modem 
microwave  communication  systems.  A  planar  array  consists  of 
a  number  of  slotted  waveguides  arranged  side-by-side,  while 
the  most  popular  way  of  feeding  the  individual  branches  is 
through  the  use  of  crossed  waveguide  couplers  with  centred- 
inclined  slots.  For  large  arrays,  the  concept  of  sub-arraying  is 
often  used  to  improve  the  usable  frequency  bandwidth,  and 
also  to  reduce  the  sensitivity  of  antenna  performance  on 
dimensional  tolerances.  Each  sub-array  is  then  fed  by  a  main 
line,  while  the  individual  main  lines  are  fed  by  a  beam 
forming  network.  Slotted  array  seeker  antennas  also  make  use 
of  sub-arraying,  where  a  monopulse  comparator  network 
effectively  feeds  the  four  quadrants  of  the  antenna  either  in- 
phase  or  out  of  phase  with  respect  to  each  other. 

Design  procedures  for  the  synthesis  of  linear  slot  arrays  were 
proposed  by  Elliott  [1-3],  and  were  later  extended  to  cover 
planar  arrays  [4].  Other  contributions  to  the  field  of  array 
design  include  those  of  [5]  and  [6].  Even  though  these 
techniques  have  been  known  and  used  for  years,  the 
inexperienced  designer  is  still  faced  with  a  number  of 
practical  problems  when  trying  to  implement  them,  especially 
in  the  case  of  larger  arrays  consisting  of  sub-arrays.  In  this 
paper,  some  helpful  hints  concerning  the  numerical 
implementation  are  provided,  and  other  practical  aspects 
which  are  not  explicitly  defined  elsewhere  in  the  literature  are 
addressed. 

The  design  procedure  involves  the  repeated  solution  of  sets  of 
nonlinear  equations  in  order  to  iteratively  compute  the  offsets 
and  lengths  of  the  radiating  slots,  as  well  as  the  inclination 
angles  and  lengths  of  coupling  slots.  The  successful  numerical 
solution  of  sets  of  nonlinear  equations  is  usually  dependent  on 
the  quality  of  the  initial  guesses  for  the  unknowns,  the  number 


of  unknowns  and  also  on  the  nature  of  the  nonlinear 
equations.  In  practice,  specifying  suitable  starting  values  may 
require  considerable  effort.  In  this  paper,  a  simple  but  highly 
effective  method,  which  is  less  dependent  on  good  initial 
guesses  and  a  priori  knowledge  about  the  eventual  geometry, 
is  proposed. 

It  is  known  that  the  solution  obtained  for  a  specific  design  is 
not  unique,  due  to  a  number  of  degrees  of  freedom.  This 
aspect  is  usually  addressed  by  fixing  some  of  the  dimensions. 
A  different  way  of  accommodating  the  redundancy  is 
proposed.  It  naturally  divides  the  nonlinear  equations  into 
smaller  groups  by  effectively  decoupling  them  into  semi¬ 
independent  sets  which  are  solved  individually  during  each 
iteration. 

There  also  seems  to  be  some  uncertainty  regarding  the  correct 
orientation  of  coupling  slots,  especially  when  the  array  has 
more  than  one  feed  line.  Unequivocal  relations  for  the  signs 
of  radiating  slot  offsets  and  coupling  slot  angles  are  thus 
provided.  Finally,  some  guidelines  on  the  specifications  for 
beam-forming  networks  required  by  multi-sub-array  systems 
are  given. 

2.  THEORETICAL  BACKGROUND 

Consider  a  planar  array  consisting  of  a  total  of  S  sub-arrays. 
The  5th  sub-array  shown  in  Fig.  1  consists  of  a  total  of  Ts 
branch  lines,  while  the  rth  branch  line  has  a  total  of  Ns  t  slots. 
The  slots  are  resonantly  spaced,  i.e.  d  =  Xg/2  and  dQ  «  X  / 4 
at  the  design  frequency,  f0.  The  inter-slot  spacing  should  be 
larger  than  the  waveguide  width,  a,  and  for  a  waveguide  filled 
with  a  dielectric  material  of  relative  permittivity  er,  the  value 
of  a  should  be  between  the  limits  X0/2^er  <,  a  <,  X0/^Ter, 

where  X0  is  the  free  space  wavelength.  The  thickness  of  the 
septum  separating  neighboring  branches  is  thus  t0  =  X/2  -  a. 
The  waveguide  height,  fc,  may  be  chosen  arbitrarily.  The  sub¬ 
array  is  fed  by  means  of  a  main  line,  which  is  connected  to  the 
different  branch  lines  via  centred-inclined  coupling  slots.  The 
coupling  slot  feeding  the  rth  branch  line  of  the  5th  sub-array 
has  an  inclination  angle  of  Qs  t  and  a  slot  length  of  l$  r  The 
angle  0s  t  is  taken  as  positive  in  a  clockwise  direction.  This 
coupling  slot  is  located  between  the  £  >;th  and  the  r+l)th 
slot  of  the  rth  branch  line.  The  term  us  =  us\  is  the  unit 

vector  in  the  direction  toward  the  shorted  end  of  the  main  line, 
so  that  us  =  ±1  for  sub-arrays  fed  from  the  bottom  and  the 
top,  respectively.  The  complete  array  has  a  total  number  of 

m  =  Y,Y,  Kt  slots- 

S=\  t= I 


1054-4887  ©2000  ACES 


28 


ACES  JOURNAL,  VOL  15,  NO.  1,  MARCH  2000 


Two  distinct  radiating  slot  indexing  conventions  are  used: 

(a)  A  local  numbering  system  with  a  triple  index  ( s ,  t,  n ) 
denoting  the  «th  slot  in  the  rth  branch  line  of  the  5th  sub- 
array.  The  slot  has  a  width  of  w,  an  offset  of 
relative  to  the  center  line  of  the  branch  line,  and  a  slot 
length  of  Lstn.  Slots  in  a  branch  line  are  numbered  from 
left  to  right,  irrespective  of  whether  the  sub-array  is  fed 
from  the  bottom  (as  in  Fig.  1)  or  from  the  top.  Branch 
lines  are  numbered  such  that  the  line  closest  to  the 
feeding  end  is  denoted  by  t  =  1  and  the  line  closest  to  the 
shorted  end  of  the  main  line  by  t  =  Ts . 

(b)  A  global  numbering  system  with  a  single  index,  i.  The 
slot  with  local  index  (5,  t9  n)  has  a  global  index  of 

S-i  Tu  /- 1 

'  =  EE  K., +  E  +  »• 

«*=  1  t-l  V  =  1 


Branch  line  Ts 
d  \ 


=>=  1 


Branch  line  t  (s,  r,  kx,)i  k?"  !  ( s ,  r, 

= _ &0”)  _  T3  M  \  =P 


*1' 


-4 


c=jCr> 


Branch  line  2  (-s',  2,  ks2) i  j 

(J»  2»  1)  4  1  c =n 

=  r0|  d\  !  j(j,  2,^+1) 

(*/j2  »  Zll ) 

Branch  line  1  t 

Tgrco"! 

\(s,  i,  *,.,+!) 

foj  » Zs,\  ) 


voltage  in  the  equivalent  network  is  Vstn.  The  term 
y5se/fw  =  yself(x^w,I5 1  n)  represents  the  normalized  self¬ 
impedance  of  this  slot.  The  function  g..  is  given  by  [2],  which 
is  dependent  on  the  coordinates  of  the  centers  of  the  rth  and 
the  yth  slots.  Let  (xi ,  zf)  the  coordinates  for  the  rth  slot  in  the 
global  coordinate  system.  In  terms  of  the  local  coordinates,  it 
is  given  by 


(xnzi)  =  (</  +  Xs°Z  >  <r  ~do~  Ns,td  +  nd)  (3) 


where  (x/( ,  z/,)  are  the  coordinates  of  the  center  point  of 
the  right  edge  of  the  rth  branch  line.  The  functions  fs  and  hstn 
may  be  obtained  from  [2]  and  [3]  by  making  the  substitutions 
21  n  =  Lstn  and  xn=x°^n.  For  n  =  lor  n  =Nsl,  the 
undefined  terms  in  (2)  such  as  Ffs'“0,  hs  t  0  and 

hstN  + ,  may  be  found  through  the  application  of  image 
theory!  The  constants  K{,  K2  and  K3  may  also  be  obtained 
from  [2]  and  [3]. 


dQ  d  d  d0 


Fig.  2  Equivalent  transmission  line  network  for  the  rth  branch 
line  of  the  5th  sub-array. 


Let  the  voltages  across  the  shunt  elements  in  Fig.  2  be  given 
by 


s,t,n 


9s.,  <,(~VKr 


(4) 


Fig.  1  Geometry  of  the  5th  sub-array. 


The  equivalent  network  for  the  rth  branch  line  of  the  5th  sub¬ 
array  is  shown  in  Fig.  2,  where  ysatn  is  the  normalized  active 
admittance  of  the  «th  slot.  The  design  procedure  involves  the 
repeated  solution  of  a  set  of  nonlinear  equations  in  order  to 
determine  the  slot  offsets  and  lengths.  The  nonlinear  equations 
are  related  to  the  so-called  design  equations,  given  by  [2] 
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The  term  qst  is  an  arbitrary  constant  equal  to  ±1,  depending 
on  whether  the  offset  of  the  slot  with  index  (s,t,ksJ)  is  to  be 
positive  or  negative.  These  voltages  are  arbitrarily  chosen  to 
be  real  valued.  Then  is  a  positive  real  number 
representing  the  magnitude  of  the  voltages  on  the  equivalent 
circuit  of  the  rth  branch  line  of  sub-array  5. 


The  series  impedances  of  the  equivalent  circuit  for  the  main 
line  in  Fig.  3  are  given  by  [4] 

z(i'°  =  (K (5) 


where  yf£’'^  is  the  normalized  input  admittance  of  the  rth 
branch  line  of  the  5th  sub-array,  given  by 


<*,o  =  y  y« 

/in  SS't'fi 


(6) 


n  -  1 


For  the  slot  with  local  index  (5,  t,  n )  and  global  index  /,  the  slot 
voltage  is  denoted  as  Vffn  or  vfot9  while  the  corresponding 
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z(s,  1)  z(s,  2)  z(s,Ts) 


d  d  d 

Fig.  3  Equivalent  network  for  the  main  line  of  the  5th  sub¬ 
array. 


The  term  is  the  coupling  coefficient  defined  by 


while  S^’0  =  Sn(Q  -  |  05J  ,  /  =  /5 ,)  is  the  scattering  matrix 
component  of  the  crossed  waveguide  couplers  formed  at  the 
junction  of  the  5th  mainline  and  the  rth  branch  line,  as  shown 
in  Fig.  4.  The  normalized  input  impedance  of  the  5th  main  line 
is  then  given  by 


(^.0)2^0 


(8) 


/=i  /=i 


2 


1 


Fig.  4  Crossed-guide  coupler  geometry  for  the  main  line  of 
the  ith  subarray  and  the  rth  stick. 


Due  to  the  fact  that  the  coupling  slots  are  also  resonantly 
spaced  (i.e.  d=  XJ2  at  the  design  frequency,  fQ),  the  current  /, 
in  Fig.  3  is  given  by 


2  a<*> 

v^(4,)+n 


where  Is  x  is  the  current  at  the  first  series  element,  while 
is  the  incident  wave  intensity.  Suppose  the  input  impedance  of 
each  sub-array  is  purely  resistive.  Then  =  r^\  where 
is  the  chosen  normalized  input  resistance  of  the  5th  sub-array 
as  seen  from  the  first  series  impedance  in  the  main  line 
equivalent  circuit.  The  relation  between  the  current  at  the 


series  impedance  representing  the  rth  coupling  slot  in  the  main 
line,  Is  t,  and  the  voltage  across  the  ks  ,th  element  in  the 
equivalent  network  for  the  rth  branch  line,  Vs  {  ^  ,  is  given  by 

0  *  S> 

V...  =  u,  .  J—  k{sJ) I  .  rim 


One  design  objective  is  to  have  all  the  slot  voltages  in  phase. 
From  equation  (1),  it  follows  that  the  argument  of  the  term 
VgtJfstn  should  be  uniform.  In  agreement  with  the  choice 
of  real-valued  equivalent  network  voltages,  this  phase 
reference  is  taken  to  be  zero.  In  a  specific  branch  line,  this  is 
accomplished  by  alternating  the  offsets  of  the  slots  on 
opposite  sides  of  the  center  line.  The  sign  of  the  term  fStitk  { 

will  be  equal  to  the  value  of  qs  r  which  implies  that  the  sign 
of  Vs  t  k  should  also  be  identical  to  the  value  of  q  r  This 
agrees  with  the  choice  in  (4).  For  resonant  coupling  slots, 
is  purely  real,  while  (9)  and  (10)  dictate  that 
should  also  be  real  valued.  Consequently,  the  orientations  of 
the  coupling  slots  are  given  by 


A  comparison  of  equations  (4),  (9)  and  (12)  then  yields  the 
simple  relation  of  [4],  given  by 


V°  =  7° 

s7t 


(13) 


where  7°  =  |  Is  , 


3.  DERIVATION  OF  THE  SETS  OF  NONLINEAR 
EQUATIONS 

Consider  the  5th  sub-array.  Select  the  ks  ;th  slot  branch  line 
t  to  serve  as  the  reference  slot.  From  equations  (1)  and  (4),  it 
is  found  that 


y  °  f  FS‘°‘  (-1  )"■*'• 

ss,tyn  _  J$yl,n  s,t,n  V  J 
a  r  TrSlot 


Note  that  the  ratio  in  ( 14)  is  only  dependant  on  the  slot  offsets 
and  lengths  in  branch  line  t  of  sub-array  5.  Equation  (14)  does 
not  provide  a  sufficient  number  of  equations  to  uniquely  solve 
for  the  unknown  slot  offsets  and  lengths  in  the  branch.  It  has 
been  suggested  that  the  slot  which  is  to  have  the  greatest  slot 
offset  be  singled  out  and  assigned  an  offset  near  the  upper 
limit  of  reliable  input  data  [4].  However,  it  is  often  difficult  to 
identify  this  slot  without  some  degree  of  a  priory  knowledge 
of  the  eventual  geometry,  especially  in  cases  of  uniform 
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excitation.  A  different  approach  is  to  specify  the  input 
admittance  of  each  branch  line  to  be  purely  real  and  with  a 
normalized  input  admittance  of  gf£’t\  Equations  (6)  and  (14) 
may  then  be  used  to  construct  the  following  set  of  2Ns  t 
nonlinear  equations  for  the  unknown  slot  offsets  and  lengths 
of  branch  line  t  in  sub-array  s : 

/  f  Kslot  (-1 

ss,t,n  _  Js,t,n  s,t,n  V  1 ) 


Re 


a  ,  p-slot 

n=l,2,..N  n*k 


Im[^J=0  n  =  \,2,..Nsf 

E  y?,t,n 


(15) 


~  Re 


W  =  1 


=  0 


where  yf  n  is  calculated  using  equation  (2). 

Comparing  the  active  admittances  of  the  ks  ,th  slot  of  the  rth 
branch  line  and  the  ks  x  th  slot  of  the  first  branch  line,  using 
equations  (1)  and  (12),  results  in 


ysaxKl  "  « 


As>{) 


(16) 


If  the  slot  offsets  and  lengths  of  the  reference  slots  are  known, 
equations  (8)  and  (16)  may  be  used  to  construct  the  following 
set  of  2  Ts  nonlinear  equations  for  the  magnitude  of  the 
inclination  angles  and  coupling  slot  lengths  of  main  line  s : 


Re 


ySjx 


I  fs 


s,t,ks 


Fslot 


V  °  I  f  I  Ksl0t 

Im  [  ]  =  0  t=\,2,..Ts 


=  0  t  =  2,3,..Ts 
(17) 


Re 


.« _  y:  (k <*■'>) 


/  =  i 


=  0 


where  is  calculated  using  (7). 


4.  EFFECTIVE  IMPLEMENTATION  OF  THE 
STARTING  VALUE  PROBLEM 

Depending  on  the  nature  of  a  set  of  nonlinear  equations,  their 
numerical  solution  may  be  highly  dependent  on  the  quality  of 
the  initial  guesses  for  the  unknowns,  as  well  as  the  number  of 
unknowns.  In  the  previous  section,  the  equations  were 
subdivided  into  smaller  sets,  but  the  designer  is  still  faced  with 
finding  good  starting  values  of  the  unknown  slot  dimensions. 

For  larger  arrays,  the  evaluation  of  the  external  mutual 
coupling  terms  g..  is  generally  the  most  time  consuming  step 
in  the  numerical  implementation  of  a  design  procedure.  For 


practical  purposes,  an  iterative  approach  for  the  design  of 
linear  arrays  was  proposed  in  [2],  where  during  each  iteration, 
these  terms  are  computed  using  the  current  set  of  slot 
dimensions  and  are  assumed  to  remain  constant  while 
calculating  new  dimensions  by  way  of  solving  of  a  set  of 
nonlinear  equations.  It  was  also  suggested  in  [4]  that  the 
external  and  internal  mutual  coupling  terms  be  set  to  zero  for 
the  first  iteration.  Due  to  the  simplified  nature  of  the  resulting 
nonlinear  equations,  this  negates  the  need  for  good  initial 
guesses  for  the  unknown  as  required  by  nonlinear  equation 
solvers.  Arbitrary  starting  values  within  the  realistic  range  of 
slot  dimensions  may  then  be  specified,  while  the  nonlinear 
equation  solver  routine  is  usually  successful  in  finding  a 
solution.  During  the  second  iteration,  the  computed  slot 
dimensions  for  the  zero  mutual  coupling  case  are  used  to  re¬ 
evaluate  the  mutual  coupling  terms,  and  the  same  dimensions 
are  employed  as  starting  values  for  the  subsequent  solution  of 
the  nonlinear  equations.  This  approach  is  entirely  effective  for 
linear  arrays.  This  is  due  to  the  fact  that  the  mutual  coupling 
between  slots  that  are  approximately  axially  aligned  is 
relatively  small,  and  therefore  the  abrupt  inclusion  of  the 
mutual  coupling  contributions  during  the  second  iteration  only 
has  a  secondary  effect,  and  the  changes  in  slot  dimensions 
never  become  excessive. 

For  planar  arrays  the  displacement  between  certain  slots  is 
also  lateral,  which  causes  the  external  mutual  coupling  to 
have  a  more  pronounced  effect  on  the  slot  dimensions.  The 
examples  in  [3]  bare  evidence  of  this.  Even  though  a  uniform 
excitation  for  all  slots  was  specified,  the  mutual  coupling 
causes  large  variations  in  slot  offsets,  while  in  the  absence  of 
mutual  coupling,  all  slots  have  the  same  lengths  and  offset 
magnitudes.  Implementation  of  the  approach  proposed  for 
linear  arrays  may  result  in  numerical  instability.  Problems 
generally  do  not  occur  when  solving  for  the  zero  mutual 
coupling  case  during  the  first  iteration.  However,  during  the 
second  iteration  where  the  zero  mutual  coupling  dimensions 
(which  may  differ  appreciably  from  the  final  results)  are  used 
to  first  calculate  g..  and  then  also  as  starting  values  for  the 
subsequent  solving  of  the  nonlinear  equations,  the  procedure 
often  fails. 

It  is  therefore  proposed  that  the  procedure  be  modified  by 
retaining  the  first  iteration  with  zero  mutual  coupling,  but 
instead  of  abruptly  introducing  the  effects  of  mutual  coupling 
during  the  second  iteration,  that  its  effects  be  included 
gradually.  This  is  achieved  by weighing  the values  of  g..  and  hstn 
with  a  weighting  factor  increasing  from  0  to  1  during  the  first 
half  of  a  chosen  total  number  of  iterations.  Consequently,  the 
first  iteration  will  be  performed  with  zero  mutual  coupling, 
while  the  full  effect  of  the  mutual  coupling  is  only  included 
during  the  second  half  of  the  iterations.  This  approach 
generously  compensates  for  a  possible  increase  in  the  total 
number  of  iterations  by  being  very  robust  and  reliable.  The 
ease  of  implementation  of  this  procedure  is  also  an  advantage 


J.C.  COETZEE,  J.  JOUBERT:  DESIGN  OF  PLANAR  SLOT  ARRAYS  REVISITED 


31 


over  other  approaches  which  merely  rely  on  good  starting  P 

values  that  are  obtained  by  for  instance  first  considering  an  ( |  a^ht\2  -  |  ^spllt  |2  )  =  0  (23) 

infinite  array  [6]. 


5.  BEAM-FORMING  NETWORK 


An  array  consisting  of  a  number  of  sub-arrays  would  naturally 
require  a  power  splitter  or  beam-forming  network  to  distribute 
the  power  from  the  antenna  input  to  the  multiple  sub-array 
main  lines.  A  complete  design  should  therefore  also  include 
the  specifications  of  the  power  splitter  network.  A  comparison 
of  the  active  admittances  of  the  k  .  th  slot  of  the  first  branch 

s,  1 

line  of  sub-array  s  and  of  the  kx  x  th  slot  of  the  first  branch  line 
of  sub-array  1  using  equations  ( 1 )  and  (12)  yields  the  following 
relative  amplitudes  of  the  incident  wave  intensities  for  the  sub¬ 
arrays: 


flW|  _  C*,,,  ^  Vi.^M) 

"  IWI  V?«K^y;xk'(r"+ 1) 


,0)| 


(18) 


The  phase  difference  between  a^s)  and  a{l)  is  either  0°  or 
180°.  The  phase  relation  is  chosen  in  accordance  with  the 
phase  properties  of  the  beam-forming  network  to  be  used. 
Together  with  (18),  the  ratio  a^/a^  is  thus  completely 
specified.  In  general,  an  array  consisting  of  S  sub-arrays  would 
require  a  splitter  with  P  =  S  +  1  ports.  For  the  sake  of  clarity, 
it  is  assumed  that  each  of  the  first  S  ports  is  connected  to  its 
corresponding  sub-array,  while  the  /Th  port  serves  as  the 
common  feed  port.  From  the  definition  of  the  scattering 
parameters,  it  follows  that 

[lS'split]  [asplit]  =  [isplit]  (19) 


[Ssplit]  is  the  P*P  scattering  parameter  matrix  of  the 
network,  specified  with  the  phase  reference  of  port  1  to  port  S 
at  the  center  of  the  first  coupling  slot  of  the  main  line  it  is 
connected  to.  Since  the  first  S  ports  of  the  splitter  network  are 
terminated  in  the  input  impedances  of  the  sub-arrays,  it 
follows  that 


split  _  p(s)  »  split 
as  1  in 


s  =  1 , 2, . .  S 


where 


-n^') 
1  in 


(0  _ 
in 


(20) 


(21) 


The  wave  emanating  from  the  Jth  port  is  fed  directly  into  the 
main  line  of  the  5th  sub-array,  while  the  antenna  should 
normally  be  matched  at  port  P.  This  implies  that 


^  split 


s  -\,2,.  .S 
[0  s=P 


(22) 


For  a  lossless  splitter  network,  the  conservation  of  power 
condition  requires  that 


Using  (22)  and  (23)  together  with  the  specification  for 
flCO/^O)  then  gives 


,  split 


split 


- en 

n(  1) 


N 


s 

E 

m= 1 


a 


(m) 


,0) 


1  -  p(m)2 


s  =  l,2,..S 

(24) 


[0  s  =  P 

C  is  an  arbitrary  common  phase  term.  The  relations  in  (24) 
specify  the  requirements  for  the  splitter  network.  The  simplest 
way  of  realizing  this  would  be  to  design  the  sub-arrays  to  be 
matched  (i.e.  r ^  =  1  and  =  0),  which  would  then  require 

that 


(q(j)/q(1))  e * 


5r  Split 

sP 


(a(m)/a(l))2 


0 


5  =  1,2,  ..S 

(25) 


s=P 


6.  DESIGN  PROCEDURE 

Array  design  requires  data  on  the  self  admittance  of  isolated 
radiating  slots  as  a  function  of  the  slot  offset  and  length,  as 
well  as  scattering  parameters  for  coupling  slots  in  terms  of  the 
inclination  angle  and  slot  length.  These  may  either  be 
obtained  through  measurement  or  calculation.  In  the  case  of 
the  latter,  the  formulations  in  [7]  and  [8]  may  be  used  to  pre 
compute  the  databases,  while  a  bivariate  spline  interpolation 
scheme  as  in  [9]  has  proved  to  be  very  effective  for  the 
calculation  of  the  terms  =  ysc,f(x°^n  >Ls  t  n)  and 

K(J-0  =  k(  |  0  ,  | ,  /  r ).  The  design  procedure  then  consists 
of  the  following: 

1 .  Specify  the  slot  voltages  vf°xn  for  a  desired  radiation 
pattern. 

2.  Specify  us  for  each  sub-array,  as  well  as  qs  l  and  ks  t 
for  each  branch  line. 

3 .  Select  realistic  values  for  the  input  conductance  of  each 

branch  line,  A  guideline  is  to  choose  atypical  slot 

offset  and  to  identify  the  corresponding  resonant 
length  for  this  offset,  L  .  A  typical  value  for  the  input 
conductance  would  then  be  gP’1*  =  Ns  l  .yMlf(x^,£.  ). 

4.  Specify  the  input  resistance  of  each  main  line,  rP\ 
Unless  the  circumstances  dictate  otherwise,  the  most 
convenient  choice  would  be  =  1 . 

5.  Select  an  appropriate  number  of  iterations  to  be 
performed, C ,  where  C  ^  10  is  an  even-valued  integer. 

6.  Select  initial  guesses  for  the  lengths  and  offsets  for  the 
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7. 

8. 

9. 

10. 
11. 


12. 

13. 

14. 


15. 

16. 


17. 

18. 
19. 


20. 


radiating  slots.  An  effective  choice  is  to  select 


off 


-*./.»  =  311(1  Ls,t,n  =  Ltyp ■ 

Select  initial  guesses  for  the  magnitude  of  the  inclination 


angles  and  slot  lengths  for  the  coupling  slots. 

Set  the  iteration  counter  c  =  1 . 

Calculate  the  global  slot  coordinates  using  equation  (3). 
Calculate  the  mutual  coupling  terms  g..  and  hs  t  n . 

For  the  first  half  of  the  total  number  of  iterations 
(c  <  C/2%  weigh  the  values  of  gJf  and  hs  t  n  with  a 
factor  of  2(c  - 1  )/C,  i.e.  scale  all  the  elements  of  g..  to 
2  ( c  ~  1 )  gjj  /  C  and  the  elements  of  hstn  to 
2{c-\)ijC. 

Set  s  =  1 . 

Set  t  =  1. 

Solve  the  set  of  2 Ns  t  nonlinear  equations  in  (15)  for  the 
unknown  slot  offsets  and  lengths  of  branch  line  t  in  sub- 
array *:  (xftn> Lstn)  where  n  =  \,2,..Nst. 
Increment  t  and  repeat  step  14  until  t  =  Ts. 

Solve  the  set  of  2TS  nonlinear  equations  in  (17)  for  the 
magnitude  of  the  inclination  angles  and  coupling  slot 
lengths  in  main  line  s :  ( I  1 1  >  4  t)  where 

t  =  l,2,..Ts. 

Increment  s  and  repeat  steps  13-16  until  s  =  S. 
Increment  c  and  repeat  steps  9-17  until  c  =  C. 

Calculate  the  ratios  of  the  wave  intensities  that  a  power 
splitter  network  needs  to  supply  to  the  different  sub¬ 
arrays  from  (18). 

Specify  the  phase  relation  at  each  of  the  first  S  ports  of 
the  power  splitter,  i.e.  whether  sub-arrays  are  to  be  fed  in 
phase  or  180°  out  of  phase  relative  to  the  excitation  of 
the  first  sub-array,  and  use  (11)  to  determine  the 
orientation  of  the  coupling  slots. 


7.  DESIGN  EXAMPLE 

As  an  example,  a  small  array  similar  to  the  one  shown  in  Fig.  5 
was  designed.  It  consists  of  four  sub-arrays,  each  having  three 
branches  and  three  slots  per  branch.  A  design  frequency  of 
f0- 9  GHz  was  chosen,  and  half-height  X-band  waveguide  of 
width  a  =  22.86  mm  and  height  6  =  5.08  mm  was  used  for  the 
branch  guides  and  feed  lines.  A  database  for  the  properties  of 
isolated  radiating  slots  and  inclined  coupling  slots  was 
computed,  using  a  slot  width  of  w  =  1.5875  mm  and  a 
waveguide  wall  thickness  of  t  =  1.27  mm.  A  uniform 
excitation  for  all  radiators  was  adopted.  The  array  was 
designed  to  have  a  normalized  conductance  of  =  1 .5  for 
all  branches,  and  a  normalized  input  resistance  of  =  1  for 
each  of  the  four  sub-arrays.  From  the  geometry,  it  follows  that 
U\  =  u2  =  ”1 ,  =  w4  =  1 ,  kl  {  =  k^  t  =  1 ,  k2  t  -  £3 ;  =  2  and 

qs  t  =  -1 .  A  total  number  of  ten  iterative  steps  were  used 
during  the  design.  The  computed  slot  offsets  and  slot  lengths 
for  the  radiating  slots  are  shown  in  Table  I,  while  the 
inclination  angles  and  lengths  of  the  coupling  slots  are  listed 
in  Table  II.  Due  to  the  symmetry,  sub-arrays  land  3  and  sub¬ 
arrays  2  and  4  are  mirror  images  of  each  other,  and  therefore 
corresponding  slots  will  have  the  same  magnitude  of  offset  or 


inclination  angle,  but  with  a  difference  in  sign. 


(2,1,1)  !=□ 

(1,1,1)  C=1  (1,1,3) 

1 _ -i  (1,1,2)  r— 1  "i 

C=D  1 - 1 

n==D 

OAD  Er?  (1A3) 
UA2)  1 - 1 

1=1  (2,3,3) 

<w>  g 

r~— )  cz=3 

1 - ,  (3,3,3) 

1=1  (4,3,3) 

1 - 1  C=1 

CTTTD  CZZZI 

CZ3 

CZZD  1 - 1 

1 - 1  1 - 1 

(3,1,1)  (=3 

(4,1,1)  1=^ 
rr-m  1 - 1 

c 

c 

Fig.  5  Geometry  of  a  4><3><3  array. 


Table  I  Slot  offsets  and  slot  lengths  of  radiating  slots. 


Slot  index  (s,t,n) 

xsl„  (mm) 

4,/,»  W 

(1,1,1),  (3,1, 3) 

t2.256 

16.758 

(1,1, 2),  (3,1,2) 

±3.015 

16.821 

(1,1,3),  (3,1,1) 

*2.397 

16.779 

(1,2,1) ,  (3,2,3) 

*2.859 

16.822 

(1,2,2) ,  (3,2,2) 

±1.958 

16.287 

(1,2, 3),  (3,2,1) 

*2.221 

16.574 

(1,3,1) ,  (3,3,3) 

*2.177 

16.842 

(1,3,2) ,  (3,3,2) 

±1.516 

16.787 

(1,3,3) ,  (3,3,1) 

*2.325 

16.830 

(2,1,1) ,  (4,1,3) 

±2.848 

16.828 

(2,1,2),  (4,1,2) 

*2.012 

16.630 

(2,1,3),  (4,1,1) 

±3.138 

16.923 

(2,2,1),  (4,2,3) 

±1.976 

16.437 

(2,2,2) ,  (4,2,2) 

*2.679 

16.610 

(2,2,3) ,  (4,2,1) 

±2.153 

16.500 

(2,3,1) ,  (4,3,3) 

±1.969 

16.888  . 

(2,3,2) ,  (4,3,2) 

*2.188 

16.756 

(2.3.3) .  (4.3. 1) 

±1.603 

16.885 

Table  II  Inclination  angles  and  slot  lengths  of  coupling  slots. 


Slot  index  (  s,  t) 

0f>,  (degrees) 

4.,  W 

O.l),  (3,D 

±23.61 

16.986 

I  d,2),  (3,2) 

*21.47 

16.970 

(1,3) ,  (3,3) 

±18.58 

16.948 

(2,1),  (4,1) 

±24.73 

16.998 

(2,2),  (4,2) 

*20.81 

16.964 

12.3) .  14.3) 

±17.88 

16.943 
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Due  to  the  symmetry,  the  incident  wave  intensities  were  found 
to  be  a^/a^  =  1  for  all  sub-arrays. 

The  effect  of  the  gradual  inclusion  of  the  mutual  coupling  is 
illustrated  in  Fig.  6,  where  the  variation  of  the  offset  for  the 
slot  with  index  (1,3,2)  as  computed  during  the  iterative  steps 
is  shown.  The  initial  offset  calculated  for  the  case  with  zero 
mutual  coupling  is  3.082  mm,  while  the  final  offset  is  only 
1.516  mm.  The  gradual  introduction  of  the  mutual  coupling 
contributions  over  the  next  five  iterative  steps  avoids  abrupt 
changes  in  the  slot  dimensions,  and  convergence  is  easily 
achieved  during  the  second  half  of  the  iterations. 


Number  of  iterations 


Fig.  6  Offset  of  slot  (1,3,2)  as  computed  during  each  itarative 
step. 

8.  CONCLUSION 

The  repeated  computation  of  the  external  mutual  coupling 
between  slots,  which  involves  numerical  integration,  is  the 
most  time-consuming  part  of  the  procedure.  Utilization  of  an 
alternative  expression  for  the  coupling  coefficient  [10]  results 
in  an  appreciable  increase  in  efficiency.  Some  numerical 
experimentation  may  be  required  in  order  to  obtain  a  design 
where  all  the  slot  offsets  are  within  the  desirable  range.  Slot 
offsets  of  less  than  w/2  should  be  avoided.  Care  should  also  be 
taken  not  to  allow  the  offsets  to  become  too  large,  especially 
when  using  waveguides  of  reduced  height.  This  may  lead  to 
problems  related  to  the  validity  of  the  equivalent  circuit 
models  for  the  slots  [11].  If  the  computed  offsets  of  a 
particular  branch  are  either  too  large  or  too  small,  the  input 
conductance  of  those  branches  should  be  decreased  or 
increased  and  the  design  repeated.  Used  in  conjunction  with 
an  analysis  procedure  [12],  the  admittance  levels  of  the  branch 
lines  may  be  chosen  so  as  to  optimize  the  off-centre  frequency 
performance  of  the  array. 

The  proposed  algorithm  is  easily  translated  into  computer 
code,  and  has  been  applied  successfully  to  the  design  of  arrays 
of  varying  sizes.  The  process  of  subdividing  the  nonlinear 
equations  into  smaller  groups  and  gradually  introducing  the 


effects  of  mutual  coupling  has  proven  to  be  very  reliable 

when  used  in  conjunction  with  reliable  nonlinear  equation 

solver  routines. 
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Note:  Code  (or  algorithm)  capability  descriptions  are 
not  acceptable,  unless  they  contain  sufficient  technical 
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Where  possible  and  appropriate,  authors  are  required  to 
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that  used  in  this  text.  A  double  column  format  similar  to 
that  used  here  is  preferred.  No  author’s  work  will  be 
turned  down  once  it  has  been  accepted  because  of  an 
inability  to  meet  the  requirements  concerning  fonts  and 
format.  Full  details  are  sent  to  the  author(s)  with  the  letter 
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